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ABSTRACT 

We report the results of an investigation of particle acceleration and electron-positron plasma 
generation at low altitude in the polar magnetic flux tubes of Rotation Powered Pulsars, when 
the stellar surface is free to emit whatever charges and currents are demanded by the force- 
free magnetosphere. We apply a new ID hybrid plasma simulation code to the dynamical 
problem, using Particle-in-Cell methods for the dynamics of the charged particles, includ- 
ing a determination of the collective electrostatic fluctuations in the plasma, combined with 
a Monte-Carlo treatment of the high energy gamma rays that mediate the formation of the 
electron-positron pairs. We assume the electric current flowing through the pair creation zone 
is fixed by the much higher inductance magnetosphere, and adopt the results of force-free 
magnetosphere models to provide the currents which must be carried by the accelerator. The 
models are spatially ID, designed to explore the physics, although of practical relevance to 
young, high voltage pulsars. 

We observe novel behavior, a) When the current density j is less than the Goldreich- 
Julian value (0 < j/ j a] < 1), space charge limited acceleration of the current carrying beam 
is mild, with the full Goldreich- Julian charge density being comprised of the charge density 
of the beam, co-existing with a cloud of electrically trapped particles with the same sign of 
charge as the beam. The voltage drops are on the order of mc 2 /e, and pair creation is absent, b) 
When the current density exceeds the Goldreich-Julian value (j/ j G! > 1), the system develops 
high voltage drops (TV or greater), causing emission of curvature gamma rays and intense 
bursts of pair creation. The bursts exhibit limit cycle behavior, with characteristic time scales 
somewhat longer than the relativistic fly-by time over distances comparable to the polar cap 
diameter (microseconds), c) In return current regions, where j/ j GJ < 0, the system develops 
similar bursts of pair creation. These discharges are similar to those encountered in previous 
calculations of pair creation when the surface has a high work function and cannot freely emit 
charge, Timokhin (2010). In cases b) and c), the intermittently generated pairs allow the sys- 
tem to simultaneously carry the magnetospherically prescribed currents and adjust the charge 
density and average electric field to force-free conditions. We also elucidate the conditions 
for pair creating beam flow to be steady (stationary with small fluctuations in the rotating 
frame), finding that such steady flows can occupy only a small fraction of the current density 
parameter space exhibited by the force-free magnetospheric model. The generic polar flow 
dynamics and pair creation is strongly time dependent. The model has an essential difference 
from almost all previous quantitative studies, in that we sought the accelerating voltage (with 
pair creation, when the voltage drops are sufficiently large; without, when they are small) as a 
function of the applied current. 

The ID results described here characterize the dependence of acceleration and pair cre- 
ation on the magnitude and sign of current. The dependence on the spatial distribution of the 
current is a multi-D problem, possibly exhibiting more chaotic behavior. We briefly outline 
possible relations of the electric field fluctuations observed in the polar flows (both with and 
without pair creation discharges) to direct emission of radio waves, as well as revive the pos- 
sible relation of the observed limit cycle behavior to microstructure in the radio emission. 
Actually modeling these effects requires the multi-D treatment, to be reported in a later paper. 
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1 INTRODUCTION 

Young pulsar wind nebulae (PWNe) show that rotation powered 
pulsars (RPPs) have dense magnetospheres, at least with regard 
to those regions that feed the plasma outflow (e.g., Bucciantini 
et al. 2011). Electron-positron pair creation in the open field line 
region that connects to the external world is the only known candi- 
date for the origin of such outflows, with acceleration and convert- 
ible gamma ray emission occurring either at low altitude (Sturrock 
1971) or in the outer magnetosphere (Cheng et al. 1986). High den- 
sity flows that can feed all the open field lines can exist only in the 
low altitude polar cap region (for a general discussion, see Arons 
2009). 

Theoretical studies of charged particle flow from the mag- 
netic polar regions of rotation powered pulsars began with the 
observation by Goldreich & Julian (1969), that an isolated mag- 
netic rotator in vacuum must have a charged magnetosphere al- 
most corotating with the star. Since RPPs are strongly gravitation- 
ally bound and cool objects (thermal scale height in any atmo- 
sphere orders of magnitude less than the stellar radius) and have 
no external source of plasma supply (so far as we know), the only 
plasma source is extraction of charged particles from the stellar 
surface, leading to a conjectured magnetosphere whose plasma is 
fully charge separated, in contrast to all other known astrophys- 
ical systems, whose plasmas are charged but quasi-neutral. Gol- 
dreich & Julian (1969) speculated that on polar field lines - those 
that extend beyond the light cylinder located at cylindrical radius 
R LC = cP/2n =a 48, 000 P km, P = rotation period in seconds - 
a charge separated outflow would form. They argued that the en- 
ergy/particle in the outflow would be no more than the gravita- 
tional escape energy GMJR, ~ O.3mc 2 (M„/1.4M )(lO km/tf,), 
M, and R, are star's mass and radius correspondingly - the parti- 
cles leave at non-relativistic energies, in spite of the fact that the 
electric potential drop across the polar field lines is equal to the full 
potential of an open rotating magnetosphere with a dipole magnetic 
field O m = y/WJc' ~ 10(P/10- IS ) [/2 P- 312 TV, W R = rotational en- 
ergy loss rate, with P = dP/dt . O m vastly exceeds the rest energy 
and gravitational energy of the particles, either electrons or protons 
(or He or other ions populating the star's crust and atmosphere). 
The super strong magnetic field suppresses free acceleration of the 
particles in the transverse electric field, whose primary ("zeroth or- 
der") consequence is corotation of the field lines with the magnetic 
field embedded in the neutron star (NS), with field line motion mea- 
sured by the E x B drift of charged particles across the magnetic 
field (which occurs even when the particles have zero Larmor gy- 
ration). The particle loss rate in the conjectured charge separated 
scenario is N R = c<t> m /e * 2 x 10 30 (M0~ 15 ) 1/2 P~ 3/2 s _1 , orders 
of magnitude less than that inferred from the injection of plasma 
into the young PWNe. The electrodynamics of the magnetosphere 
differs drastically, depending on whether the particle loss rate falls 
short of or exceeds N R . For the young PWNe, the particle injection 
rate exceeds N R (by a lot). In that case, the magnetosphere's basic 
state should be one in which E ■ B = 0, with no parallel acceleration 
sufficient to generate convertible gamma rays occurring under the 
pulsar's rotational control. 

The discovery of gamma ray pulsars in the 1970s, and their 
proliferation into a population with more than 100 such stars in the 
most recent published Fermi pulsar catalog (Abdo et al. 2010), has 
shown that parallel acceleration to GeV gamma ray emitting en- 
ergies (indeed, multi-hundred GeV, in the Crab pulsar, Aliu et al. 
2011) must occur somewhere, with energy efficiency exceeding a 
few tenths of a percent, as measured by L y /W R , L y = gamma-ray 



luminosity. If the acceleration is limited by radiation reaction, as is 
true in many models, L y is a good proxy for the energy put into par- 
allel acceleration. L y /W R can approach as much as 50% at smaller 
spin down luminosities. Just how some fraction of the total poten- 
tial drop gets released in acceleration along B, gamma ray emission 
and pair creation has been mysterious since the beginning of pulsar 
research, made relevant to the real world by the gamma ray discov- 
eries. Since the polar cap source is the only one capable of feeding 
the whole (open) magnetosphere, its understanding remains of cen- 
tral interest to modeling pulsar magnetospheres, even though the 
spectral and beaming characteristics of the pulsed gamma rays are 
better modeled by accelerators in the outer magnetosphere. 

Free particle outflow from the NS surface is a common as- 
sumption in most of the current pulsar models (see §2). The polar 
cap accelerator problem has been studied under that assumption 
before (e.g. Michel 1974; Fawley et al. 1977; Mestel et al. 1985; 
Shibata 1997; Beloborodov 2008). Michel (1974) and Fawley et al. 
(1977) obtained solutions for the non-neutral space charge limited 
charged particle flow for the current density almost equal to the 
GJ current density. All these models assume strictly steady flow 
in the co-rotating frame, on all time scales. In these models, the 
charge density of the current carrying beam supplies almost all of 
the charge density needed to short out the parallel component of 
the electric field, while leaving a residuum E\\ sufficient to acceler- 
ate the beam - relativistic energies in a temporally steady flow are 
found if the current density j\\ = -(B/P) cos^ + small corrections = 
jai\ X = Z (A*> the pulsar inclination angle. Mestel et al. (1985) 
showed that the velocity of the beam is monotonically increasing 
with altitude to relativistic speeds only if the current density is 
larger than j OI . If the current density is smaller that j GS the tempo- 
rally steady velocity of the beam (assumed to have no momentum 
dispersion) oscillates spatially, i.e. particles accelerate and decel- 
erate to a complete halt as they move outwards into the magneto- 
sphere. Beloborodov (2008) rediscovered Mestel et al.'s solution 
and suggested that in the region of the polar cap where y'y < j GJ 
particles will be not accelerated up to high energies as the beam ve- 
locity oscillates, but along the magnetic field lines with j 01 > 
or j\\/j G j < particle acceleration will be efficient and will lead to 
pair formation. The quantitative model which we describe in this 
paper lends some support to Beloborodov's speculations, although 
it does not agree with them in detail. 

In this paper we describe our study of the physics of the po- 
lar cap accelerator in the space charge limited flow regime start- 
ing from first principles - assuming free particle outflow from the 
surface of a NS we compute the electric field, particle accelera- 
tion, gamma-ray emission, propagation and pair creation simul- 
taneously. It extends the study of current flow and pair cascades 
in neutron star magnetospheres using the theoretical formulation 
and self-consistent numerical techniques introduced in Timokhin 
(2010). 

The plan of the paper is as follows. In §2, we review the prop- 
erties of the current flow imposed by the magnetosphere, in the 
force-free model, pointing out that the current density is the main 
parameter which regulates the efficiency of particle acceleration. In 
§3 we review the properties of stationary solutions for the charge 
separated space charge limited flow problem. In §4 we briefly de- 
scribe our numerical model. In §5 we describe the results of nu- 
merical modeling for the case of sub-GJ current density, the regime 
when particle acceleration is inefficient and no pair creation is pos- 
sible. In §6 we consider flow regimes with efficient pair creation. In 
§6.4 we pay special attention to the stationary flow regime, which 
up to now was assumed in (most) works on pulsar polar cap accel- 
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erators, and describe why it has limited relevance to the force-free 
model of the pulsar magnetosphere. We discuss the implications of 
our results for the physics of rotation powered pulsars in §7 and 
summarize our conclusions in §8. 



2 CURRENT DENSITY IN THE POLAR CAP 

Nebular observations of plasma supply by RPPs suggest the open 
field line regions are "MHD-like", i.e., having E\\ = except in 
special zones (such as the polar cap), which are, in effect, bound- 
ary layers. There is essentially no observational information on the 
properties of the closed magnetosphere, thus the simplest hypothe- 
sis is to follow Goldreich & Julian (1969) and assume the magneto- 
sphere is (almost everywhere) filled with plasma that shorts out Ey. 
The magnetosphere has open and closed magnetic field line zones. 
In the open zone plasma flows away into the pulsar wind; currents 
and their associated electromagnetic inertia keep the magnetic field 
open, so sustaining the flow. The plasma flows all the way from 
the base of the open field line zone in the polar cap of the pulsar 
where it is either extracted from the NS surface or generated by 
electron-positron cascades (or both). The distribution of the current 
density across the open field line zone, and therefore, across the po- 
lar cap, is determined by the global magnetospheric structure. Sta- 
bility of the pulsar mean profiles and sharpness of the peaks in the 
spectra of gamma-ray pulsars strongly suggest that on scales com- 
parable to the light cylinder, the magnetosphere is stationary in the 
frame co-rotating with the neutron star with smooth and continuous 
plasma outflow. However, the stationary co-rotating magnetosphere 
hypothesis demands stationarity only in a statistical sense, fast lo- 
cal fluctuations which average to a stationary state (and even more 
broadly, global variations) can be included within this picture, so 
long as they do not smear the beaming profiles. 

The polar cap acceleration and possible pair cascade zone - 
the main place where electron - positron plasma that feeds the wind 
can be produced - is much smaller than the characteristic scale R LC 
of the magnetosphere. Therefore its inductance is negligible com- 
pared with that of the magnetosphere, and the polar current flow 
within the polar cap region must have an average equal to that set 
by the magnetosphere's global structure. 

In this paper we solve a local problem of how the polar cap 
cascade zone adjusts itself to the current density required by the 
magnetosphere. On the dynamical timescales typical for the cas- 
cade zone (microseconds) the magnetospherically required current 
density is stationary because it could change only on magneto- 
spheric time scales (tens of milliseconds up to several seconds). 
The idea that the acceleration zone has a magnetospherically deter- 
mined current appears in the electrodynamics through the magnetic 
induction equation 



at 

with 



= -4*ry|| + c(V x B\ * -An(J\\ ~ jm), 



jm — ^ ^magnetosphere)|| 



(i) 



(2) 



being the current that sustains the twist to the field lines. In this 
paper we neglect the induced variations in the magnetic field that 
accompany variable E\\ - because of the very strong background 
magnetic field (which in the region of interest has V x B = 0), 
these have little dynamical effect on the acceleration and cascade 
dynamics, see Appendix A for derivation of eq. (1). 



We study the behavior of the cascade zone under different cur- 
rent loads, sampling the range of possibilities illustrated in Fig. 1. 
The model has an essential difference from previous studies, in that 
we seek the accelerating voltage (with pair creation when the volt- 
age drops are sufficiently large, without when they are small) as a 
function of the applied current. Previous work has almost entirely 
focused on the opposite direction, seeking the current that emerges 
from the accelerator when the voltage is fixed, either by the ge- 
ometry of by the poisoning of the accelerator by pair creation. In 
addition, we allow the system to be fully time dependent. These 
generalizations lead to qualitatively different results from what has 
appeared before. Our model is one-dimensional, with spatial axis 
along magnetic field lines; from here on we drop subscript || from 
all quantities. 

The characteristic charge density, the Goldreich- Julian charge 
density, 



n b 

2nc 



sets the characteristic current density 



Jgi — t]Gi£ ■ 



(3) 



(4) 



Following Fawley et al. (1977) and Arons & Scharlemann (1979), 
we assume also that if the star lacks an atmosphere, the work func- 
tion for charged particles to leave the NS surface is small enough 
that any number of charged particles can be extracted from the NS 
surface until the extracting electric field is screened. In the classi- 
cal pulsar regime, the theory of charged particle binding to the crust 
suggests such free emission to be likely (Medin & Lai 2010). More 
relevantly, X-ray observations of heated polar caps (Bogdanov et al. 
2007, and references therein) suggest these stars have atmospheres 
overlying the solid and ocean components of the crust, which guar- 
antees free emission of charges. The opposite case, with complete 
suppression of particle emission from the surface, was studied in 
Timokhin (2010), a realization of the scenario conjectured by Rud- 
erman & Sutherland (1975) 

As we will show in the next sections there are three quali- 
tatively different plasma flow types in the polar cap of pulsar de- 
pending on the ratio of the current density imposed by the magneto- 
sphere to the GJ current density: (i) j m has the same sign and its ab- 
solute value is smaller that the GJ current density, < j m / j Q1 < 1, 
hereafter sub-GJ current density; (ii) j m has the same sign and its 
absolute value is larger that the GJ current density, j m /j GJ > L 
hereafter super-GJ current density; (iii) j m has the opposite sign to 
the GJ current density, j m / j G] < 0, hereafter anti-GJ current den- 
sity. 

The advent of quantitative solutions for the structure of the 
force-free model (e.g. Contopoulos et al. 1999; Timokhin 2006; 
Spitkovsky 2006; Kalapotharakos & Contopoulos 2009; Bai & 
Spitkovsky 2010) has provided, for the first time, a theory of the 
current flow expected as a function of pulsar parameters (S, P, x, 
when the magnetic field is a star-centered dipole). Earlier modeling 
of polar cap, slot gap and outer gap accelerators adopted the expec- 
tation that the current density is on the order of j al , and expressed 
the hope that the accelerator and pair creation physics do not sen- 
sitively depend on the precise value and spatial distribution of j. 
The results reported here show that the magnitude and sign of the 
current flow do lead to drastic differences in the open field line ac- 
celerator's behavior in the three regimes (i) - (iii), even though the 
order of magnitude of the current is as expected. Fig. 1, showing 
J = 3 1 jai, reveals that all three flow regimes occur in the force-free 
magnetosphere model. While \J\ always has numerical values on 
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Figure 1. Field-aligned current density at the polar cap of the force-free rotator, with j\\ measured in units of the Goldreich-Julian current density ja = 
-£1 ■ BI2nc. The black circle is the rim of the polar cap - the footprints of the field lines that pass outside the light cylinder fall with that circle. The distributed 
current is shown. The current sheet component coincides with the polar cap boundary. This plot was made by Xuening Bai using results of force-free 
magnetosphere simulations presented in Bai & Spitkovsky (2010). 



the order of unity, it can be negative (return current) as well as ly- 
ing in the separate regimes < J < 1 and j > 1. We show these 
separate regimes have different dynamical behavior and different 
implications for pair creation. 

We also show that once the constraints of the steady flow mod- 
els for space charge limited flow are relaxed, the small departures 
of the GJ charge density from the simple estimate -f26cos^/27rc 
created by geometric and general relativistic considerations (Arons 
& Scharlemann 1979; Muslimov & Tsygan 1992) that play an es- 
sential role in the steady flow models 1 have little significance when 
the flow is fully time dependent. Since we consider only the polar 
cap region, with altitudes not exceeding the width of the polar flux 
tube r pc = R«iR,/R LC = 0.145.P~ 1/2 km <k R, = 10 km, spatial 
variation of B is mostly unimportant. If we do not say so explic- 
itly otherwise, throughout the paper we assume that the GJ charge 
density is constant, independent of the distance along B. 



3 TIME-STATIONARY SPACE CHARGE LIMITED 
FLOW 

Copious pair creation occurs when there is a sufficiently large ac- 
celerating potential difference along B when pairs are absent. Such 
potential drops, typically > 10 12 Volts (Sturrock 1971), readily ex- 
ist in the absence of current flow, that is, in a vacuum. Under pulsar 
conditions, if current does flow, pairs appear when the current flow 
co-exists with TV potential drops 2 , that is the current is a relativistic 
beam. Therefore, the starting place is the properties of time station- 
ary, space charge limited, charge separated flow. In this section we 
give an overview of these properties while detailed derivation of 

1 If the GJ charge density were uniform and the beam is everywhere rela- 
tivistic and time stationary, the unique model is no acceleration at all (Tade- 
maru 1973; Fawley et al. 1977), a severe contradiction. 

2 TV potential drops are required if curvature radiation from charges ac- 
celerated along a locally dipole magnetic field is the emission mechanism 
for the gamma rays that convert to pairs. The required potential drops are 
smaller if the emission mechanism is inverse Compton scattering (either 
resonant or non-resonant) of softer photons from the surface (e.g. Hib- 
schman & Arons 2001, and references therein). 




x [X D ] 

Figure 2. Phase space trajectories (4-velocity p vs. distance normalized to 
^d. oi) for stationary space charge limited flow for current densities j m f j as = 
0.1, 0.25, 0.5, 0.75, 0.9, 0.95 

the equations used in this section is given in Appendix B together 
with some useful asymptotics. These review and extend a variety 
of results already in the literature. 

For definiteness, we consider pulsars with an acute angle x 
between SI and B ("acute" pulsars). These objects have t] G! < at 
the polar cap, and require electron emission to supply r] GJ . 

Let us consider an electron beam starting with zero velocity at 
x = and let the current density j m imposed by the magnetosphere 
(eq. (2)) be a fraction £ of the GJ charge density 

Jm = £/oj • (5) 

In our case the GJ charge density is negative. In stationary flow the 
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Figure 3. Phase space trajectories for stationary space charge limited flow 
for current densities j m /jci = 0.9,0.95,0.99, 1, 1.1, 1.5. Note that in con- 
trast to Fig. 2 the vertical axis on this plot is logarithmic. 
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Figure 4. Charge density of stationary space charge limited flow normalized 
to the GJ charge density r/ GS as a function of distance x normalized to A D for 
current densities j m / j G1 = 0.5, 1.0, 1.5. 



current density is constant in both space and time and is equal to 
the imposed current density j = j m , thus dE/dt = in eq. (1). The 
stationary electric field £ s is then given by Gauss's law which in 
the frame corotating with the NS takes the form (e.g. Fawley et al. 
1977, and references therein) 



dE s 
dx 



■ 4n(r] - ?/ G j) . 



(6) 



The magnitude of the electric field increases with distance if the 
magnitude of the charge density rj is larger than the GJ charge den- 
sity and decreases otherwise. 

The charge density at any given point of the flow is 



n = j/v ■■ 



Jn 



P 



(7) 



where v is the flow velocity and p = yv/c is the 4-velocity = mo- 
mentum in units of inc. p of a charge separated stationary beam is 
given by the solution of the equation (B10) 



dp] 2 ~ P 2 + 1 



ds 



2^(1+^- VFTT) 

pi \ i 



(8) 



where the distance s is measured in units of the Debye length /1 D . c 
of a cold electron plasma with GJ number density 

v-l/2 



4TTT] Gl e 



(9) 



B\2 is the pulsar magnetic field in 10 G, P is the pulsar period in 
seconds. 

The numerical solutions of (8) for different imposed current 
densities are shown in Figs. 2, 3 for different values of these 
results confirm the work of Shibata (1997). For < j m / j G! < 1 the 
steady flow oscillates spatially, with particle momenta oscillating 
in the interval [0, p max ], with 



Pu 



(10) 



dp /ds = at p = p max and so is the RHS of eq. (8). The value of 
Pmax and the spatial period of oscillations s increase with increas- 
ing see Figs. 2, 3. For £ > 1 acceleration is monotonic with p 
increasing to infinity, with asymptotic behavior 



P 



V2.5- 



f- 1 



(11) 



See Fig. 3 for£ = 1, 1.1, 1.5. 

The reason for such behavior is as follows. The flow starts 
with zero initial velocity at x = where the electric field is zero. 
When the particle velocity is small the imposed current density is 
produced by high particle density moving slowly. In such places 
the absolute value of the beam charge density (cf. eq. 7) is larger 
that that of the GJ charge density, \rj\ > \rj GJ \ - see Fig. 4 where 
we plot the ratio rj/r] GS for flows with £ = 0.5, 1, 1.5. The charge 
density is negative and according to eq. (6) the electric field in this 
region is decreasing towards more negative values, thus accelerat- 
ing electrons. If the imposed current density exceeds the GJ current 
density, f > 1, the absolute value of the beam charge density - 
whose maximum value is \j m /c\ = £\r] GI \ - never becomes smaller 
than I^qjI, hence, dE s /dx < and acceleration continues up to in- 
finity - the electric field and potential are monotonic. For f < 1, 
on the other hand, particles' velocities increase to values such that 
the imposed current density can be sustained by particle number 
density smaller than the GJ number density. Then \t]\ < \r] as \ (see 
the line for j m = 0.5j G! in Fig. 4), dE s /dx > and the accelerat- 
ing electric field weakens, changes sign and decelerates particles - 
in cold flow, the particles decelerate to zero velocity, and the cycle 
repeats. 

The model of space charge limited flow outlined here pro- 
vides a physical framework for the expected particle energetics. It is 
based on the approximation of one cold fluid and an assumption of 
complete stationarity of the flow. It can be extended to a two-fluid 
model to account for presence of positrons and pair creation (as 
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in Arons 1983). However kinetic effects such as particle trapping 
can not be included in a cold fluid approximation, although certain 
aspects can be modeled if momentum dispersion ("pressure") is in- 
cluded, with an assumed equation of state. A kinetic model incor- 
porates momentum dispersion in the collisionless medium without 
having to make arbitrary assumptions about the equation of state. 
As we show in the following sections, particle trapping and pair 
creation profoundly affect the plasma dynamics, with momentum 
dispersion being essential to the dynamics behind the simultaneous 
adjustment of the charge density to the condition of low voltage 
drop along B, modeled as E ■ B = in the force- free global model, 
and the adjustment of the field aligned current j to the magneto- 
spherically imposed j m . 

The study of plasma kinetics in a general regime - without 
relying on stationarity or perturbation theory - is possible only by 
means of numerical simulations. In following sections we describe 
our study of plasma, both fully non-neutral and quasi-neutral when 
pair cascades form, with the help of a self-consistent hybrid numer- 
ical model incorporating both charged particles and photons code. 



4 NUMERICAL SETUP 

We use the same one-dimensional hybrid Particle-In-Cell/Monte 
Carlo hybrid code described in Timokhin (2010) modified for the 
space charge limited flow regime. Below we briefly describe the 
main equations, notations and numerical algorithms; a detailed de- 
scription can be found in Timokhin (2010), §§2, 3. 

We solve the evolutionary equation for the electric field E 

dE(x, t) 



dt 



-4n(j(x,t)- j m ) 



(12) 



-47r(?7 s tart " ^Gj) 



(13) 



(14) 



where j(x, t) is the actual current density and j m is the current den- 
sity imposed on the cascade zone by the magnetosphere. This equa- 
tion is the Ampere's law, eq. (1). We are solving an initial value 
problem, thus an initial distribution of the electric field E(x, t = 0) 
must be supplied. At the start of the simulation we construct the 
initial distribution of the electric field by solving the Gauss equa- 
tion for the electric potential <p assuming some initial charge density 
distribution ^ stan 

dx 2 

E = -** 

dx 

We proceed with the time integration of eq. (12) using a charge 
conserving scheme (e.g. Villasenor & Buneman 1992; Birdsall & 
Langdon 1985), so the Gauss equation is satisfied at each succes- 
sive time step up to machine precision. The GJ charge density en- 
ters in the equation (13) for the initial configuration of the electric 
field; this information is then "carried on" in time by eq. (12). 

To model the space charge limited flow at every time step we 
inject electrons and protons just outside the numerical domain used 
for electric field calculation and let the system pull the necessary 
amount of particles into that domain. We do not set E(x = 0,0 
to zero as a boundary condition but rather allow the plasma in the 
system to enforce this condition as part of the simulated physics. 
A detailed description of our algorithm for reproducing the space 
charge limited flow condition at the NS surface is given in Ap- 
pendix C. When pair creation cascades occur, we take into account 
only curvature radiation as the gamma-ray emission mechanism; 
pairs are created by single photon absorption in the strong mag- 
netic field (e.g. Erber 1966). 



We performed many numerical experiments starting from dif- 
ferent initial conditions: (i) computational domain filled by plasma 
with charge density equal, less, and higher than the GJ charge den- 
sity as well as starting with vacuum; (ii) different initial potential 
drop over the domain; (iii) different length of the computation do- 
main. In all cases without exception after initial relaxation on the 
time scale of the order of the flyby time of the domain the system 
settled down to a configuration which depends only on the imposed 
current density j m . 



5 LOW ENERGY CHARGE SEPARATED, SUB-GJ 
FLOW: < jJj G} < 1 

As it turns out, there is no pair formation for < j m / j a < 1 and 
the only characteristic spatial scale of the flow is the Debye length 
/1 D G ,. In this section we will discuss the properties of such flow us- 
ing simulations in domains with the length L up to several hundreds 
of 

-^d,gj» which is much less than the width of the polar cap. Using 
such a small domain we well resolve the characteristic spatial scale; 
increasing the domain length does not change the results. 

According to the stationary solution from §3 there is no rel- 
ativistic particle acceleration when < j m /jos < 1- The flow is 
spatially oscillatory and the maximum momentum of particles p m . ix 
is by far too small for emission of pair producing photons. It is 
unlikely that such oscillatory cold flow can exist - near the stagna- 
tion points it may be a subject to "instability" with characteristics 
of the wave breaking to which nonlinear waves in cold fluids with 
velocity stagnation are subject, which could destroy the spatial os- 
cillation and create an effective resistor, across which a substantial 
fraction of the perpendicular (to B) voltage drop might appear. The 
available voltage drop across B is huge and in principle it might 
happen that the system ends up in a state with highly oscillatory 
electric field and bursts of pair formation as predicted in the two- 
fluid plasma model of Levinson et al. (2005). On the other hand, 
finite amplitude electrostatic waves containing similar stagnation 
points show wave breaking, with alteration of cold to hot flow with 
momentum dispersion no more than comparable to the flow veloci- 
ties found in the originally constructed cold oscillation (Akhiezer & 
Polovin 1956; Tajima & Dawson 1979). That would create a warm 
but non-relativistic or mildly relativistic charge- separated outflow. 

We find that in the sub-GJ regime, the non-neutral flow is in- 
deed low-energy, with particle energies orders of magnitude below 
the energy/particle required for pair production. However, the final 
state differs drastically from the oscillatory flow from the stationary 
solution shown in Figs. 2, 3. Even if at the beginning of the sim- 
ulations particles follow trajectories of the stationary solution, the 
standing non-linear wave structure breaks quickly - particles at the 
velocity zeros of the cold flow go both up and down. 

A good example of this inherent instability is shown in Fig. 5. 
We start from a vacuum configuration - when there are no particles 
in the domain - and let the system evolve. In Fig. 5 we show snap- 
shots of the phase space portraits of the flow with j m = 0.5 j C i, the 
distance from the NS, normalized to A D G] , is along x-axis, particle 
momenta, normalized to m e c, are along y-axis. The whole domain 
has the length L = 50A D GJ and we show only a part of it here. Time t 
is measured in flyby time of the whole domain L/c. With the dashed 
red line we show phase space trajectories of particles from station- 
ary solution. Particles coming from the surface at first follow the 
trajectories of the stationary solution. However, after coming to the 
first stagnation points some of the particles are turned back and the 
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Figure 5. Development of "trapped particle" flow when the space charge limited flow which starts into vacuum. Phase space portrait of particles are shown for 
16 moments of time indicated in small boxes on the top of each plot. The current density j m = 0.5 j GI . The distance x is measured in units of the Debye length. 
Red dashed lines show the analytical solution for stationary flow. Particle momenta p_ are normalized to m e c. The total length of the computational domain 
L = 50/l D . ai, only part of it (x < 50A D G j) is shown here. Time is measured in flyby time (Lie) of the whole domain. Snapshots in the same row have the same 
time interval between them, but these time intervals are different for different rows, they increase towards the bottom row. 



flow start to randomize . After several tens of plasma periods A Di G1 /c 
the flow reaches its final configuration. 

Examples of final configurations for space charge limited flow 
with different current densities are shown in Fig. 6, where we plot 
phase space portrait of particles in the whole computation domain. 
The flow has two components: a warm beam of particles with high- 
est momentum which produce the required current density and a 
cloud of charged particles circulating in the domain - these com- 
pose an electrically trapped, "thermal" component. In the cloud 
component there is roughly equal number of particles moving in 
opposite directions, these particles do not contribute to the current 
but contribute r] GS - j m /c to the charge density keeping T] totii \ equal 
to the GJ charge density. The distinction between these components 
is not absolute as some particles from the beam go into the cloud 
and vise versa, although the fraction of mixing particles is small. 
Some of the particles in the cloud have very low momenta and so 
they can adjust to any given charge density. Hence, in sub-GJ space 
charge limited flow, < j m / j GJ < 1, the electric field is not sen- 
sitive to variation of the GJ charge density 3 - in contrast to the 



3 We performed simulations for j m / j G1 < 1 with variable GJ charge den- 
sity, and, as expected saw the electric field to be just as screened as in the 
uniform GJ density case. 



large importance of variation of the GJ charge density for relativis- 
tic acceleration of space charge limited cold beams in the polar cap 
cascade models of Arons & Scharlemann (1979) and Muslimov & 
Tsygan(1992). 

Plasma flow in the sub-GJ regime can be described as a 
beam of mildly relativistic particles propagating through a cloud of 
trapped particles with near-thermal distribution. In Fig. 7 we plot 
particle distribution functions. The beam component is visible as 
a bump on the distribution function at the high momentum side. 
The cloud component has a near-thermal (Maxwell-Juttner) mo- 
mentum distribution (at least in its low-energy part) dn/dp oc p - 
such quasi-thermalization is a common consequence of the phase 
mixing between the particles and fields built into the wave breaking 
process. 

When particles leave the NS they are non-relativistic and their 
charge density rj = j m /v = £j G! /v is larger than the GJ charge den- 
sity and so they form a charge sheet near the surface generating ac- 
celerating electric field, just as in the idealized stationary case (see 
Fig. 4). When particles reach the velocity such that \j m \/v < \r] G! \ 
the electric field derivative dE/ds changes sign and after some dis- 
tance the electric field can start decelerating particles. At the point 
where E = particles reach their maximum momentum. Above the 
gap particles from the cloud component add additional charge and 
so the change density there is equal to rj G1 and the electric field is 
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Figure 6. Phase space portraits of well developed space charge limited flows for 6 different current densities: j m / j os = 0.1, 0.25, 0.5, 0.75, 0.9, 0.95. Current 
density j m is indicated in the left upper part on each plot. Distance is measured in units of A D , particle momenta p- are normalized to m e c. Note that the lengths 
of the computation domain differ between these plots. 




Figure 7. Particle distribution functions p (dn/dp) for well developed space charge limited flows from Fig. 6. Distribution functions of particles with positive 
momenta (moving toward the magnetosphere) are shown by solid lines, distribution functions of particles with negative momenta (moving toward the NS) by 
dashed lines. 



screened. The length of this gap is the order of .s m;lx = s /2 (half the 
spatial period of cold flow oscillations) and so the maximum mo- 
mentum particles gain in this gap is comparable to p max (both s msK 
and p miix depend on j m ), In Fig. 10 we plot momenta of particles 
in the beam component as functions of £ = j m /j GS superimposed 
of the theoretical dependence of p max given by eq. (10); the agree- 
ment is pretty good. The current density in the final configuration is 
close to j m throughout the domain, with small fluctuations around 
this value Sj : <k j m . 

In Fig. 8 we plot electric field in the calculation domain at the 



same moments of time as the phase space portraits shown in Fig. 6. 
The electric field at any given point fluctuates, but the relative fluc- 
tuation at the beginning and at the end of the domain (accelerating 
and decelerating regions, see below) are much smaller than those in 
the center of the domain (the region of the charge cloud). The elec- 
tric field in the gap near the NS launches the beam component. The 
electric field at the other end of the computation domain supports 
the cloud components by reversing the momenta of most of the par- 
ticles in the cloud moving away from the NS, sending them back. 
This electric field is not strong enough to reflect most of the beam 
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Figure 8. Snapshots of electric field E for well developed space charge limited flows taken at the same moment as the phase space portraits from Fig. 6. E is 
normalized to nr] G] A D _ a ■ 
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Figure 9. Power spectra of the electric field I k = \Et\ 2 for well developed space charge limited flows from Fig. 6. k is normalized to the corresponding 1 /A D , GJ . 



particles. The mixing between the beam and cloud components is 
the strongest here. This electric field appears self-consistently when 
the flow reaches its steady configuration 4 , because it is needed to 
sustain the cloud component necessary to match both the charge 
and the current density in the domain. In our small scale ID simu- 
lations we impose the current density on the domain, which finally 
gives rise to the electric field at the right end of the domain. In real- 
ity it is the magnetosphere which sets the current density by twist- 
ing the magnetic field lines and generating electric field reversing 



4 The formation of the cloud component is not linked to the appearance of 
the electric field at the end of the domain, see e.g. Fig. 5 



some of the particles. This second region with unscreened electric 
field at the magnetosphere end of the domain in our model may be 
a "compressed" version of some parts of the outer magnetosphere. 
For example, on field lines that pass through the null surface where 
SI ■ B = 0, the plasma cloud and beam, composed of only one sign 
of charge, cannot freely enter the outer magnetosphere (Scharle- 
mann et al. 1978) in the absence of other sources of electric field 
in other parts of the magnetosphere (Goldreich & Julian's "hanging 
charge clouds"), offering the possibility of opening a vacuum gap 
within the otherwise force-free structure. On there other hand, on 
polar field lines that never cross the null surface - most of them, 
in the aligned rotator - the charge separated beam and cloud can 
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pair creation attenuation with the distance we set the magnetic field 
strength to zero starting at distance x B from the NS - the magnetic 
field is given by 



B(x). 



Figure 10. Momenta of the current currying particle beam. The widths of 
the peaks at Fig. 7 as the function of £ = j m / j G1 are shown by the blue bars. 
The solid red line is the maximum achieved momentum for the stationary 
SCLF solution. 



extend outwards "forever", in principle. Multi-dimensional parti- 
cle simulations, analogous to those of Spitkovsky & Arons (2002), 
are required to see if indeed this speculation is true (as well as de- 
termine how this essentially ID model might fit together with the 
other, more "lively" aspects of the polar flow outlined in §6). 

In Fig. 9 we show power spectra of the fluctuating electric 
field in the central parts of the calculation domain, outside of the 
acceleration zones. l k = \E k \ 2 , where E k is the spatial Fourier am- 
plitude of the electric field and the wave vector k is normalized to 
the /Iq'cj. These spectra can be fit with the power law l k oc kr a with 
a between 2 and 3 for all current densities. 

Thus, along magnetic field lines where the current density is 
sub-GJ there is no pair formation, so long as strong electric fields 
in neighboring, more active flow zones do not leak into the cloud. 
In an aligned rotator, for example, no pairs are forming above most 
of the polar cap area. Plasma flowing along such magnetic field 
lines is mildly relativistic, consists of particles of only one sign 
(electrons) and its density is low, equal to the GJ number density. 
However, observational evidence for ongoing pair formation in pul- 
sars is very strong and, hence, there must be regions in the pulsar 
magnetosphere where electron-positron plasma is created. 



6 PLASMA FLOW WITH PAIR FORMATION - 
SUPER-GJ AND RETURN CURRENT REGIONS 

Let us now consider what happens along magnetic field lines where 
the current density has either (i) the opposite sign to the GJ current 
density, // j GJ < 0, anti-GJ flow, or (ii) the same sign and absolute 
value larger than the GJ current density, jj j a > 1 , super-GJ flow. 
Case (i) includes the regions of return current, including the current 
sheet, while case (ii) is relevant for most of the magnetic field lines 
in a nearly-orthogonal rotator. As we show in this section, in both of 
these cases a strong accelerating electric field is generated. When 
the resulting potential drop is sufficient to accelerate particles up 
to the energies such that they can emit photons that convert to pairs 
within the accelerator, the plasma flow will be highly-nonstationary 
with intermittent pair creation. 

In a real pulsar the magnetic field strength falls with the dis- 
tance ( oc r~ 3 for a dipole field) and pair creation is possible only 
at sufficiently low altitude (if gamma ray interaction with thermal 
low energy photons is neglected). In order to imitate the effect of 



So, if x < x B 
0, if x > x B . 



(15) 



If the charge density were formed from steadily flowing rela- 
tivistic beams, it would vary with the height as rj(x) oc B(x). The GJ 
charge density changes in a slightly different way, rj(x) oc B(x)f(x) 
due to inertial frame dragging (Muslimov & Tsygan 1992; Beskin 
1990) or/and field line curvature (Scharlemann et al. 1978). If one 
neglects the latter effects, scaling oc B(x) can be incorporated into 
the spatial (for eq. (13)) or temporal (for eq. (12)) coordinates. The 
the electrodynamics of the cascade zone can be modeled in a ID 
problem with constant GJ charge density; the only effect of the GJ 
charge density variation will be in changing the spatial (and tem- 
poral) scales. The same scaled ID model can also be used to study 
the effects on the electrodynamics of the cascade zone of deviation 
of the GJ charge density scaling from being oc B(x) by considering 
a problem where rj GJ depends on the distance; the variation of the 
GJ charge density will be given by f(x) as the dependence on B(x) 
is already incorporated in the model. First, in §§6.1, 6.2, we con- 
sider the case when the GJ charge density is constant, i.e. this case 
corresponds to a model where variation of r] G i/B is neglected. Then 
in §6.3 we address the influence of the GJ charge density variation 
on the physics of the cascade zone. 

While simulations of space charge limited flow with sub-GJ 
current densities described in §5 can be considered as directly re- 
lated to more complete pulsar models - the distance over which 
wave breaking and trapping control the flow is much smaller then 
the width of the polar cap, and the ID approximation is well mo- 
tivated - the pair creation models presented in this section can be 
considered only as illustrative of the physics but not fully applica- 
ble to pulsars. The spatial scales over which pair creating photons 
are absorbed, even when they are emitted at very low altitude (e.g. 
the attenuation length of the magnetic field) are much larger that 
the width of the polar cap, making transverse structure essential for 
modeling the pulsar environment - such effects are necessary for 
a full evaluation of the pair yield, since much of the pair creation 
occurs in regions beyond the acceleration zone. Nevertheless, from 
our ID models we obtain insight into the basic cascades physics in 
a regime never previously modeled. Multi-dimensional models will 
be considered elsewhere. 

Within the context of the ID model, the domain length L, 
and the value of x B (the height above which the magnetic field 
is too weak to support pair creation) are the parameters having 
the largest departure from what would appear in multi-D context. 
These lengths in our simulations are much smaller that in real pul- 
sars, but that choice allowed us to construct models utilizing rea- 
sonable amounts of CPU time and to explore a broad parameter 
space. L and x B were chosen in such a way that the minimum size 
of the accelerating region necessary to start pair creation is at least 
2-3 times less than x B . 

We performed numerous simulations for different initial con- 
ditions and physical parameters in order to study qualitative be- 
havior of cascades. We performed simulations with different values 
of the numerical parameters (spatial resolution, time steps, parti- 
cle injection rate, number of particles per cell) in order to check 
the numerical convergence. In all physical cases presented in this 
section plasma flow is quasi-periodic, and this behavior does not 
depend on initial conditions - after a short relaxation time, com- 
parable to the flyby time of a relativistic particle through the com- 
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putational domain, the system settles down to a limit cycle sort of 
behavior. We describe here a particular set of simulations which 
is representative for all other models. In §§6.1, 6.2 we use simu- 
lations with the following physical parameters: the length of the 
domain L = 2.4 x 10 4 cm, the potential drop in vacuum across 
the domain AV = 10 14 Volts 5 , the radius of curvature of the mag- 
netic field lines with respect to a photon orbit p c = 10 6 cm ~ the 
stellar radius (small compared to the pure star centered dipole value 
^R,R LC ~ 10 7 8 P 1/2 cm but easily attained in offset dipole geometry 
(Arons 1998; Harding & Muslimov 2011; Arons, in preparation); 
magnetic field strength Bo = 10 12 G. The distance x B marking the 
transition to the outer magnetosphere with small magnetic field is 
set to x B = 0.7L. The particular simulations described in these two 
sections differ only in the imposed current density j m . 

We illustrate the flow dynamics with series of snapshots for 
different physical quantities shown in Figs. 11-16, 20-22. In those 
figures each column shows detailed information about physical 
conditions in the computation domain at a given moment of time: 
the number densities of electrons and positrons n+ (shown as charge 
densities of electrons and positrons r/ ± , n+ = \q ± \), total charge den- 
sity rj, current density j, the accelerating electric field E, phase 
portraits of electrons, positrons, and pair producing photons, and 
(in Figs. 14-16) - protons. Particles with positive values of the 4- 
velocity p are those which move away from the NS (toward higher 
altitude), particles with negative p move toward the NS. These plots 
are similar to ones in Timokhin (2010), with the only difference that 
now we use semi-logarithmic scale for particle momenta (linear for 
-5 < p < 5 and logarithmic everywhere else) on phase space por- 
traits that show dynamics of high and low energy particles on the 
same plot. 

The number density, charge density, and the current density 
are normalized to the corresponding GJ values: rj+ and 77 are nor- 
malized to \t]%\, j to \r/° s \c, where 77°, is the GJ charge density at 
the NS surface (the distinction between t/ g] and 77^ will be impor- 
tant in §§6.3, 6.4). The electric field is normalized to E = \rf a] \nE. 
The distance x on those plots is normalized to L, much larger than 
the Debye length A D GI . The time t is normalized to the relativis- 
tic flyby time of the computational domain L/c = 0.8 psec for the 
parameters chosen. The time is counted from the start of a partic- 
ular simulation, so only time intervals between the snapshots have 
physical meaning. 

In all cases plasma flow and pair formation has limit cycle 
behavior. In each case we illustrate this behavior by 3 series of 
snapshots taken within one typical cycle. These 3 series show 3 
phases of plasma flow: cascade ignition (Figs. 11, 14, 20), develop- 
ment of the cascade (Figs. 12, 15, 21), and filling the domain with 
dense pair plasma (Figs. 13, 16, 22). On each figure time intervals 
between snapshots are equal, but these intervals are different on 
different figures. 



5 Such vacuum potential drops over the domains of this size are realistic 
in young, high magnetospheric voltage pulsars, O m > 10 15 V. This choice 
allows studying pair formation over distances small compared to the polar 
flux time diameter. In more common pulsars with <D m < 10 13 V, pair for- 
mation happens over (much) larger distances. However, as we mentioned 
before, our simulations represent a toy model addressing the general be- 
havior of the polar cap acceleration zone and our choice of parameters is 
motivated by convenience of simulations, rather than by attempt to model a 
real pulsar. 



6.1 Flow with jJj G3 < 

In this section we describe cascade development for two cases 
when the imposed current density is anti-GJ (i.e. has the opposite 
sign to the GJ current density), for j m = -0.5 j as and j m = -1.5 j C i- 
To support such current electrons on average must move toward 
smaller x, down toward the NS, positrons toward larger x, up into 
the magnetosphere. Electrons could be freely emitted from the NS 
surface, but because the imposed current causes the growing elec- 
tric to point away from the surface, electrons at the surface are ac- 
celerated downwards and none are extracted from the star 6 . This 
situation resembles the case of the Ruderman & Sutherland (1975) 
cascades studied in Timokhin (2010), in the sense that all the par- 
ticles are produced during the burst of pair formation. Hence, the 
physics of pair cascades with anti-GJ imposed current density dis- 
cussed in this section is also applicable to Ruderman & Sutherland 
cascades as well. 

Particles leave the domain and some time after the burst of pair 
formation the domain becomes depleted of electrons. This process 
starts at the right end of the domain - at the "magnetosphere" end - 
as electrons are moving down toward the NS. In the region depleted 
of electrons, the positrons support the current density j < j m with 
charge density 77 > tj gs . This gives rise to the electric field which ac- 
celerates positrons toward the magnetosphere - see the phase por- 
traits of positrons in Figs. 11,14. The electric field in this region 
at any given point is growing with time and the size of the region 
with unscreened electric field is getting bigger as the remaining 
electrons are moving toward the NS. The electric field grows lin- 
early with the distance because positrons are relativistic and so their 
charge density remains constant (see plots for E in Figs. 11,14). 

Positrons in the region with the unscreened electric field are 
accelerated up to higher and higher energies and begin emitting 
pair producing gamma-rays. Positrons remaining from the previ- 
ous burst of pair formation are the particles which ignite the next 
burst. As positrons are moving up (and so do the first pair pro- 
ducing photons), the first pairs are produced at the largest distance 
from the NS where pair formation is still possible, in our case near 
x = x B = 0.7L. The injected pairs are picked up by a very strong 
electric field and are accelerated to high energies in less time that 
the first generation positrons. They emit pair producing capable 
gamma-rays, but now both secondary electrons and positrons are 
emitting pairs. As electrons and positrons move in opposite direc- 
tions so do the gamma-rays - see snapshots for gamma-rays phase 
space at t > 8.100 in Fig. 11 and t > 6.390 in Fig. 14. 

Secondary electrons and positrons are moving in opposite di- 
rections and the plasma gets polarized - see plots for 77+ and 77 
at t = 8.260 in Fig. 11 and at t = 6.430 in Fig. 14. When 
their number density become comparable to the GJ number den- 
sity these particles start screening the electric field (see plots for 
E at t = 8.260,8.270 in Figs. 11, 12 and at t = 6.430,6.440 in 
Figs. 14, 15). The screening starts in the region where the first pairs 
were injected because pairs have been injected here for the longest 
time and their number density is larger. 

The rate at which particles left from the previous burst of pair 
formation are leaving the domain depends on the imposed current 
density. The average bulk motion of electrons for j m = -0.5 j a] is 



6 Due to numerical noise the electric field fluctuates and sometimes elec- 
trons "from the NS surface" enter the domain; however, the number of such 
electrons is well below the fluctuation of electron number density due to 
numerical noise. 
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Figure 11. Ignition of pair formation in anti-GJ flow with y' ra = -0.5 j Gs . Several physical quantities are shown as functions of the distance x from the NS; 
x is normalized to the domain size L. Plots in each column (for the same time t) are aligned - they share the same values of x. The following quantities are 
plotted: I" row: rj ± - charge density of electrons (negative values, blue line) and positrons (positive values, red line); r\ ± is normalized to the absolute value of 
the Goldreich-Julian charge density |^ GJ |. 2 row: the total charge density r\ normalized to the absolute value of the Goldreich-Julian charge density |^ G j|- 3'"' 
row: current density j normalized to the absolute value of the Goldreich-Julian current density | j al \ = \rj G! c\. 4 rd row: accelerating electric field E normalized 
to the "vacuum" electric field Eq a \i] G! \7rL. 5 row: phase space portrait of positrons (horizontal axis - positron position x, vertical axis - positron momentum 
p + normalized to m e c). The vertical axis is logarithmic except for the region around zero momentum (-5 < p + < 5), where the scale is linear. 6* row: phase 
space portrait of electrons (horizontal axis - electron position x, vertical axis - electron momentum p- normalized to m e c). I' 1 ' row: phase space portrait of 
pair-producing photons (horizontal axis - photon position x, vertical axis - photon momentum p y normalized to m e c). 



less than 0.5c, while for j m = -1.5 j GI it is close to c 7 , and the size of 
the electron-depleted region grows slower in the former case. Sec- 
ond generation electrons, which mark the upper boundary of the 
gap, move with ultrarelativistic velocity toward the NS. Because of 
this in the case of j m = -0.5 j GJ the gap with the accelerating elec- 
tric field quickly disappears - at t = 8.430 in Fig. 12 secondary 
particles have already caught up with the particles from the previ- 



7 Note that the time intervals between the snapshots in Fig. 1 1 is two times 
larger than that in Fig. 14 



ous burst of pair formation and the electric field is screened every- 
where where pair formation is possible. In the case of j m = —1.5 j G1 
the gap moves toward the NS approximately retaining its size - see 
Figs. 15, 16. This behavior is generic. Namely, for imposed current 
densities with less-than-GJ absolute values, |y m | < |y GJ |, the average 
bulk motion of particles from the previous burst of pair formation 
is non-relativistic; that leads to quick closure of the accelerating 
gap. If the absolute value of the imposed current density is larger- 
than-GJ, \j m \ > \j GS \, the average bulk motion of particles from the 
previous burst of pair formation is relativistic and the gap propa- 
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Figure 12. Screening of the electric field in anti-GJ flow with j m = -0.5 j a ■ The same quantities are plotted as in Fig. 1 1 . 



gates in the direction of this bulk motion. The same behavior was 
observed also in the case of Ruderman-Sutherland cascades studied 
inTimokhin(2010). 

In the case of j m = -0.5 j al the gap with the accelerating elec- 
tric field closes and particles in the region where pair formation 
is possible are not accelerated anymore. The pair formation con- 
tinues only because of particles accelerated before the gap disap- 
pears. There are comparable numbers of particles accelerated in 
both directions, but electrons moving toward the NS propagate in 
the region where pair formation is possible and so after screening 
of the electric field most of the pairs are created with the initial 
momenta directed toward the NS (see phase portraits for pair pro- 
ducing gamma-rays for t > 8.180 in Figs. 11, 12, 13 and particle 
distribution functions in Fig. 17). In the case of j m = -1.5 j GJ the 
gap propagates down to the NS surface. While the gap is moving, 
positrons left from the previous bursts of pair formation, which are 
still present at the NS side of the domain, enter the gap, are picked 
up by the electric field, and emit pair producing gamma-rays. In 



this case the pair production is sustained not only by electrons ac- 
celerated in the initial screening event but also by positrons contin- 
uously accelerated in the moving gap toward the magnetosphere. 
Phase portraits of gamma-rays in Figs. 14,15,16 clearly show the 
presence of gamma-rays moving towards the magnetosphere dur- 
ing all stages of cascade development (see also the plot for particle 
distribution functions in Fig. 18). Gamma-rays with positive mo- 
menta occupy a larger and large fraction of the space as the gap 
propagates toward the NS. 

Regions far from the NS surface experience charge starvation 
first. The pair formation is possible only in the regions with strong 
magnetic field and until the starved region extends into the strong 
field domain the electric field remain unscreened. After the start of 
pair formation, the dense plasma propagates also in the direction of 
the magnetosphere screening the accelerating electric field there. 
In our simulations the region with x > 0.7L represents the rest of 
the magnetosphere, and in snapshots with t = 8.430 - 8.990, in 
Figs. 12, 13 and t = 6.560 - 6.930 in Figs. 15, 16, one can clearly 
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Figure 13. Filling of computation domain with dense pair plasma in anti-GJ flow with j m = -0.5 The same quantities are plotted as in Fig. 11. 



see how the dense pair plasma fills the domain with x > 0.7L. 
Filling of that domain with plasma is forced by a small electric field 
induced in the dense pair plasma by the imposed current. The rate 
at which these regions are filled with the dense plasma depends on 
the imposed current - in the case of j m = -0.5 j GI the average bulk 
motion of the pair plasma is about 0.5c, while for j m = -1.5 j G! it is 
relativistic. 

The accelerating electric field is very strong in the domain 
where pairs can not be produced until electrons generated in the dis- 
charge arrive, and positrons (either primary or secondary ones) en- 
tering this domain before electron arrival are accelerated up to high 
energies. In our model pairs are created only by single photon ab- 
sorption in strong magnetic field, what automatically restricts pair 
formation to the polar cap region. However electron-positron cre- 
ation via y-y absorption can be possible at much higher distances 
from the polar cap, similar to the outer gap scenario (Cheng et al. 
1986). In our model between the bursts of polar cap pair formation 
in the magnetospheric region above the low altitude pair produc- 



ing zone, at x > x B , particles can be accelerated up to radiation- 
reaction limited energies and emit very high energy gamma-rays. 
In a more realistic model, when y-y absorption is taken into ac- 
count, this should give rise to cascades in the outer magnetosphere 
resembling to some extend the outer gap scenario, except that the 
accelerating zone would not be limited to the region around the 
null surface, where the GJ charge density changes sign. Due to in- 
termittency of plasma flow in the return current region both polar 
cap and the outer gaps like cascades might exist along the same 
magnetic field lines. However, without higher dimensional simula- 
tions it is not clear how the interaction between two such cascade 
zones would play out. 

Screening of the accelerating electric field during each burst 
of pair formation gives rise to a superluminal electrostatic wave. In 
Fig. 19 we show screening of the electric field in the gap for the 
flow with j m = —0.5j O! , for flow with j m = — 1.5 \j GS the situation 
is very similar. Plots in Fig. 19 show in more detail than in Fig. 12 
how the accelerating electric field is being screened by newly cre- 
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Figure 14. Ignition of pair formation in anti-GJ flow with j m = -1.5j' GJ . The same quantities are plotted as in Fig. 1 1 with addition of the phase space portraits 
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ated pair plasma about / = 8.350. Screening of the electric field 
starts in the middle of the blob of newly created plasma and spreads 
to its edges. This spreading occurs in the form of an electrostatic 
wave which phase velocity is larger that c. This process is almost 
identical to what was observed for cascades described in Timokhin 
(2010) where more detailed discussion of the screening process can 
be found in §4.2 (cf. our Fig. 19 with Fig. 5 in Timokhin (2010)). 

As in the simulations described in Timokhin (2010), second 
generation particle momentum distributions are very broad, extend- 
ing towards non-relativistic energies. At least one of the reasons for 
momentum broadening of the particle spectra is trapping of parti- 
cles in the strong electrostatic wave excited at the beginning of each 
burst of pair formation. The dynamics of these discharges is very 



similar to what is seen in simulations of the Ruderman & Suther- 
land cascade; more details on momentum spreading can be found in 
§4.3 of Timokhin (2010). The appearance of the low energy compo- 
nent in the pair plasma is visible in the electron and positron phase 
portraits - initially at distances where pairs are injected, there are 
no particles in the phase space with low momenta, but later parti- 
cles fill in the low momentum regions (see plots at t > 8.350 in 
Figs. 12, 13 and for t > 6.500 in Figs. 15, 16). In Figs. 17, 18 we 
plot the particle momentum distribution p (dn/dp) for three dif- 
ferent moments of time. In the upper panel we plot the momen- 
tum distribution of particles moving toward the magnetosphere, p 
is positive; in the lower panel - the momentum distribution of par- 
ticles moving toward the NS, p is negative. These distributions are 
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Figure 15. Screening of the electric field in anti-GJ flow with j m = —1.5 ja- The same quantities are plotted as in Fig. 14. 



averages for regions where the electric field at the magnetosphere 
end is screened. On these plots one can see that the low energy 
plasma component is present at all stages of cascade development. 

There is another important effect in cascades along field lines 
carrying return current - positive ions (in our model protons) can 
be pulled from the NS. Both electrons and protons can be extracted 
from the NS surface, and in our simulations we allow injection 
of both species. The imposed current density requires positively 
charged particles to move toward the magnetosphere, i.e. protons 
could be pulled from the NS. However, for less-than-GJ current 
densities Qj m \ < |j GJ |) there are always electrons and positrons near 
the NS surface left from the previous burst of pair formation. The 
gap's accelerating electric field disappears before it reaches the NS, 
while pair plasma appears after every burst of pair formation keep- 



ing the electric field near NS screened. In our simulations of cas- 
cades with less-than-GJ current densities we did not see injection of 
protons. In the case of larger-than-GJ current density (]j m \ > | / GJ |) 
protons are indeed pulled from the NS, at least at some stage of 
cascade development, and this cannot be attributed only to numer- 
ical effects. In this case we see some protons pulled from the NS at 
any time, even if pair plasma is present near the NS. The number 
density of these protons never exceeds ~ 0.2« GJ , at least 10 times 
less than the number density of electrons and positrons; protons 
do not play a substantial role in the electrodynamics of discharges. 
The number density of protons pulled from the NS while electrons 
and positrons are still present near the NS surface depends on the 
numerical resolution and, hence, our simulations are inconclusive 
about whether this effect is real. However, there is a stage in cas- 
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cade development when the region near the NS gets depleted of pair 
plasma; at those times proton extraction from the NS is certainly a 
real effect. The number of particles left from the previous burst of 
pair formation near the NS is not enough to support the imposed 
current density and screen the electric field until the new portion 
of the pair plasma arrives. A gap with the electric field appears at 
the NS surface - see snapshots at t = 6.730, 6.830 in Fig. 16. This 
gap is best visible on phase portraits of positrons as a second region 
(the one close to the NS) depleted of those particles. The electric 
field in the gap starts accelerating protons from the NS surface, as 
is clearly visible in phase space portrait for protons in Fig. 16. In 
our simulations this second gap does not become large enough to 
accelerate electrons to energies sufficient to start another burst of 
pair formation near the NS surface, but we cannot exclude this pos- 



sibility in all cases. Finally, the main gap reaches the star and pulls 
out protons as well - see snapshot at t = 6.930 in Fig. 16. 

6.2 Blow with j m /j G] > 1 

In this section we describe cascade development for the case when 
the imposed current density is super-GJ. As an example of such 
flow we consider the case with j m = l.5j GS . To support this cur- 
rent density electrons must move up, towards the magnetosphere, 
and positrons down, toward the NS. There is a continuous source 
of electrons at the surface of the NS, x = 0. Positrons appear dur- 
ing the bursts of pair formation and there is no source of positrons 
during the quiet phase, between successive bursts of pair forma- 
tion. To keep the electric field screened there must be both elec- 
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Figure 17. Momentum distributions of particles at three moments of time for cascade with anti-GJ current density j m = -0.5 j al . Positron distributions are 
shown by solid blue lines, electron distributions - by dashed red lines, distribution of pair producing photons - by dotted black lines. Plots in the top row show 
distributions for particles moving toward the magnetosphere (p > 0), plots in the bottom row - distributions for particles moving toward the NS (p < 0). Each 
column corresponds to the same moment of time shown above the plots. Plots in each columns are aligned and share the same values of \p\. The following 
spacial regions were used for plotting the average distribution functions: x e [0, 0.7] for t = 8.350, x e [0, 0.75] for t = 8.5 10 and x e [0, 0.85] for t = 8.710. 




Figure 18. Momentum distributions of particles at three moments of time for cascade with anti-GJ current density j m = -1.5 j cs . Notations are the same as 
in Fig. 17. The following spacial regions were used for plotting the average distribution functions: x e [0,0.75] for t = 6.560, x e [0,0.88] for t = 6.730 and 
x e [0, 1] for t = 6.930. 



trons and positrons present at each point - electrons with larger- 
than-GJ number density moving up and positrons moving down. 
Positrons compensate for the larger-than-GJ number density of 
electrons keeping the total charge density equal to the GJ charge 
density. 

Positrons are leaving the domain and as there is no source of 
positrons, some time after the burst of pair formation the domain 



becomes depleted of positrons. This depletion starts at the upper 
end of the domain - at the "magnetosphere" - as positrons are mov- 
ing toward the NS. In the region depleted of positrons electrons pro- 
duce current density \j\ < \j m \ and the charge density |^| > |^ G ,|, be- 
cause their number density is set at regions closer to the NS where 
positrons are still present. This gives rise to the accelerating electric 
field which grows as the remaining positrons are moving toward the 
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Figure 19. Screening of the electric field and formation of superluminal electrostatic wave for cascade with j m = -0.5/ G j. There are six snapshots for the 
electric field E, the total change density and the charge density of electrons (negative values, blue line) and positrons (positive values, red line) r/ ± . All 
quantities are plotted as functions of distance x for the part of the calculation domain with intense pair formation. Snapshots are taken at equally separated 
time intervals. Plots in each column are aligned and share the same values of x. The same normalizations for physical quantities are used as in Fig. 11. The 
three thin red vertical lines in each plot mark fiducial points moving with the speed of light toward the magnetosphere. 



NS, see Fig. 20. For time shots at t = 9.380 - 9.500 on the plots 
for the charge density 77 and the current density j there are small 
jumps in both rj and j at the point where the number density of 
positrons drops to zero. The electron number density remains the 
same, but positrons are leaving, enlarging the region depleted of 
positive charges. 

Electrons are accelerated up to higher and higher energies in 
this gap and start emitting pair producing gamma-rays. In this case, 
electrons extracted from the NS surface, with the number density 
of the order of \j m \/e, are the particles which ignite the next burst of 
pair creation. As electrons are moving up (as do the first pair pro- 
ducing photons), the first pairs are produced at the farthest distance 
where pairs formation starts to be possible, near x = x B = 0.7L in 
our example. The injected pairs are picked up by the very strong 
electric field and are accelerated to high energies in less time than 
the primary electrons. Soon they start emitting pair producing capa- 
ble gamma-rays which decay into pairs, see snapshots for t > 9.420 
in Fig. 20. Electrons and positrons are moving in opposite direc- 



tions, the pair plasma becomes polarized and starts screening the 
electric field, see plots for r]+, 77, and E at t = 9.500,9.520 in 
Figs. 20,21. 

The screening of the electric field proceeds similarly to that 
in the case of anti-GJ flow with j m = —l.5j G! described in §6.1. 
The electric field is being screened first near x = x B , this creates a 
finite size gap with accelerating electric field which moves toward 
the NS. The bulk motion of positrons left from the previous burst 
of pair formation is relativistic and the newly created pairs are rel- 
ativistic too, so the boundaries of the regions where electric field 
is screened are moving in the same direction (toward the NS) with 
the same speed. The region in between with non-screened field - 
the gap - is therefore also moving, retaining its size. Electrons are 
continuously extracted from the NS to sustain the imposed current 
density. They enter the gap, are accelerated and emit pair producing 
gamma-rays. Pairs are produced by (i) primary electrons extracted 
from the NS surface accelerated by the moving gap and (ii) sec- 
ondary positrons moving toward the NS accelerated in the initial 
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Figure 20. Ignition of pair formation in super-GJ flow with j m = 1.5 y'cj- The same quantities are plotted as in Fig. 11. 



screening event. The beam of primary electrons accelerated in the 
gap is clearly seen on plots for particle momentum distribution, 
Fig. 23, as a spike in the electron distribution function for positive 
values of p (the upper plots in Fig. 23). In contrast to the case with 
j m = —1.5 jai, where the number of particles which were continu- 
ously accelerated in the gap was small and most of the pair were 
produced with momenta directed toward the NS, here more pairs 
are injected with momenta directed toward the magnetosphere - 
see Fig. 23 and plots for particle number density r]+ at t > 9.640 in 
Figs. 21, 22. As in the case of anti-GJ flow secondary particles have 
broad momentum spectra - the low energy plasma components is 
present at at all stages of cascade development, Fig. 23. 

The pair creation ends when the down flow component of the 
pairs reaches the stellar surface - the gap closes and the accelerat- 
ing electric field is poisoned throughout the strong magnetic field 
region, while the beam component that emits the primary pair cre- 
ating photons and its daughter pairs moves up at the speed of light, 
into the magnetosphere, as can be seen in the e + and e~ phase por- 



traits with advancing time in Figs. 21, 22. Positrons' bulk motion 
is toward the NS and some time after the plasma burst leaves - the 
wake of the departing burst being the source of positrons - the up- 
per regions become deplete of positrons and the cycle starts again. 

Both in the case of super-GJ flow with constant GJ charge den- 
sity described here and in anti-GJ flow described in §6.1, the cas- 
cade starts at the upper boundary of the strong field region when 
the upper regions become depleted of positrons or electrons cor- 
respondingly, as the latter move down. The rate at which these re- 
gions become charge depleted depends on the imposed current den- 
sity, the higher the absolute value of j m the faster is the bulk motion 
of the charge component moving toward the NS. For \ j m \ compa- 
rable or larger than \j G1 \ this bulk velocity is close to c. Thus the 
limit cycle period is larger than the relativistic fly-by time over the 
length of the strong field region. In our simulations, that region has 
size 0.7L = 168 meters, with fly by time 0.56/Jsec. However, in a 
real pulsar, the optical depth for pair creation by curvature photons 
emitted by the beam exceeds unity all the way out to heights com- 
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Figure 21. Screening of the electric field in super-GJ flow with j m = 1 .5j os . The same quantities are plotted as in Fig. 1 1 . 



parable to the stellar radius. Then the fly-by and recurrence times 
might be as long ~ 30 ^/sec, with the intermittent plasma starved re- 
gions taking the form of long filaments. Assessing these structures 
and accounting for competition between neighboring filamentary 
gaps is an intrinsically 2D issue, and will be treated elsewhere. 



Because the cascade starts far from the NS, at the largest dis- 
tance where pair formation becomes possible, in our simulations 
the accelerating electric field is unscreened between bursts of pair 
formation above the pair formation zone. In this regard in the case 
of constant GJ charge density super-GJ flow is similar to the anti- 
GJ flow. Namely another cascade zone in the outer magnetosphere 
might coexist with the polar cap accelerating zone. However, as 
we will show in the next section variation of the GJ charge density 
with the distance can significantly change the behavior of cascades 
in the super-GJ flow. 



6.3 Effects of spatially varying Goldreich- Julian charge 
density 

In this section we address the effect of the mismatch between the 
charge density of the space charge limited beam component and 
the GJ charge density on the cascade development. In the widely 
adopted steady flow model of polar cap cascades (Arons & Scharle- 
mann 1979; Muslimov & Tsygan 1992) this mismatch is the reason 
for the appearance of the accelerating electric field. 

We performed simulations with different scaling of the GJ 
charge density, both in shape (linear, power-law, exponential, de- 
creasing/increasing with the distance) as well as in amplitude for all 
the cases considered in the previous sections. For the space charge 
limited flow with < j m / j GJ < 1 and j m / j G! < we found no 
qualitative difference in cascade behavior due to variation of the 
GJ charge density. Namely, for anti-GJ flows, j m /j G s < 0, the re- 
gions far from the neutron star became charge-depleted first and the 
cascade started at the largest distance where pair formation is possi- 
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ble, with the acceleration zone propagating toward the NS. Sub-GJ 
flows, with < j m l j GJ < 1, — if the imposed current density was 
less that the local value of the GJ current density rj GI (x)c every- 
where - remained low energetic and had two-component: a moder- 
ately relativistic beam propagating through a cloud of trapped par- 
ticles. The main properties of space charge limited flow for these 
cases described in §§5, 6.1 are not affected by variations of the GJ 
charge density with the distance. 

The flow with super-GJ current density j m /j os > 1, however, 
can be strongly affected by variations of the GJ charge density - 
we observed different flow behavior when the absolute value of the 
r] GS /B increased with the distance from the NS. In that case the gap 
can appear close to the NS surface, in contrast to all other cases 
when it appeared at large distances. If such gaps can generate pairs, 
the regions above such gaps, far from the NS, remain filled with 
dense pair plasma at all times. 

Below we describe in detail an example of such flow using 
simulations with exaggerated charge density contrast in order to 



clearly demonstrate the above mentioned effects. We assume that 
the GJ charge density changes with the distance x as 



T](X) : 



1 + 



•III 



(16) 



where 77^ is the GJ charge density at the surface of the NS and a is 
a positive number. 

Very close to the surface (x <k R,) all the sources of inhomo- 
geneity of t] G! / B can be approximated in this manner. For example, 
if variations of the GJ charge density are because of the field line 
curvature (Arons & Scharlemann 1979) the parameter a in eq. (16) 
for dipolar magnetic field would be given by the following expres- 
sion 

3 0, cos^.sin^ L 

a AS ~ , (17) 

4 cos^- - (3/2) 6, cos <f>„ sin^- R, 

where 6, and (/>, are colatitude and azimuth of the field line at the 
NS surface, cf. eq. (12) in Arons & Scharlemann (1979). Note that 
for favorably curved magnetic field lines cos cf>, < and so a AS > 0. 
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Figure 23. Momentum distributions of particles at three moments of time for cascade with super-GJ current density y' ra = 1.5 j aI . Notations are the same as 
in Fig. 17. The following spatial regions were used for plotting the average distribution functions: x e [0,0.795] for 1 = 9.640, x e [0, 1] for t = 9.880 and 
.re [0,1] fort = 10.100. 



If GJ charge density variations are due to inertial frame dragging 
(Muslimov & Tsygan 1992) then 

L 

a MT x 3k s cosx— , (18) 

Kg = 2GI/RIC 2 ~ 0.15, / = NS's moment of inertia, cf. eq. (32) in 
Muslimov & Tsygan (1992). 

Physical parameters of our simulations are chosen in such a 
way that the length of the accelerating gap is smaller that x B . This 
setup differs from those described in §§6.1, 6.2 in that the vacuum 
potential drop over the the domain is larger - AK = 9 x 10 14 Volts 8 . 
The parameter in the expression (16) for the GJ charge density is 
a = 0.7, so that the GJ charge density varies linearly from — I^7gj I 
to -Ul^l throughout the domain 9 . The imposed current density 
is j m = 2 jgj = 2rj%c, so that everywhere in the domain the im- 
posed current density is larger than the local value of the GJ cur- 
rent density, j m > ri 0J (x) c. As the variation of the GJ charge density 
with the distance (after accounting for decrease in B, see the second 
paragraph of the §6) are not larger than 15% (Hibschman & Axons 
2001), our setup can be considered as a toy model for cascades in 
young pulsars with j m > 1.15 . All other parameters of the model 
are the same as those of the models in §§6.1, 6.2. 

In contrast to all other cases considered in previous sections 
the dense pair plasma is always present at large distances from the 
NS - the intermittent temporal gaps separating periods when the 
domain is full of plasma disappear. So, we illustrate our example 
with two series of snapshots in Figs. 24, 25. 

When the whole domain is filled with pair plasma the cur- 
rent density is constant, j = j m . The number density of electrons 



see footnote on p. 5 
9 This amount of variation is huge compared to what actually occurs - 
for example, if the variation is due to inertial frame dragging, for realistic 
parameters a ~ 0.01, if L ~ r pc . 



is larger than \ij GJ \/e and positrons (which are left from the previ- 
ous burst of pair formation) are necessary to keep the charge den- 
sity equal to the local value of the GJ charge density. The absolute 
value of the GJ charge density is smaller closer to the NS surface, 
Tj al is negative, and more positrons are necessary to screen the elec- 
tric field there than at larger altitudes. Positrons, on average, move 
toward the NS, and, as they slam on the NS surface and leave the 
domain, positrons from higher altitudes should come closer to the 
NS in order to keep the electric field screened. At some altitude the 
required flux of positrons towards the NS becomes larger than the 
positron flux that can be extracted from an adjacent higher altitude 
region without making that region charge starved - the resulting 
charge starvation causes a starvation electric field to appear close 
to the NS (see snapshot at t = 9.528 in Fig. 24, the gap is best visi- 
ble in the phase portraits on electrons and positrons). In the region 
above the gap there are still enough positrons to screen the electric 
field - less positrons are needed here as the absolute value of the GJ 
charge density is higher there and so the mismatch in charge density 
ly'm/c - r](x)\ is smaller. At lower altitudes, after all positrons which 
had been there before the gap appeared leave the domain, the influx 
of positrons from above is unable to keep the electric field screened 
and so the gap extends up to the NS surface (snapshot at t = 9.636 
in Fig. 24). 

The upper end of the gap extends toward higher altitudes as 
positrons move toward the NS, depleting higher altitude regions of 
positive charge. While the gap develops, until the ignition of the 
new cascade, the current density remains close to the imposed cur- 
rent density - electrons are being extracted from the NS surface and 
provide the required current density (see plots for j in Fig. 24). The 
current density is super-GJ, j m > j a , there are too few positrons 
in the gap and \rj\ > \t] gj \ what gives rise to a strong accelerating 
electric field in the gap (see plot for E at t = 9.636 in Fig. 24). 
Both electrons extracted from the NS surface and positrons flowing 
from higher altitudes are being accelerated in the gap and as the 
gap widens their maximum energies get higher and they become 
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Figure 24. Ignition of pair formation in the flow with linearly varying GJ charge density and the imposed current density j m = 2j GJ . The same quantities are 
plotted as in Fig. 24 



the primary particles igniting the new discharge cycle (snapshot at 
( = 9.744 in Fig. 24). 

Filling of the gap with pair plasma proceeds quickly, as the 
pairs are injected throughout the whole length of the gap (time shots 
at t = 9.780, 9.816 in Fig. 25). This happens because both electrons 
and positrons, moving in opposite directions, are emitting pair cre- 
ation capable photons, and their number densities are comparable. 
The gap is filled with plasma before the particles created in the pre- 
vious burst of pair formation leave the higher altitude region, and 
so, in contrast to the cases described in § §6. 1 , 6.2 the higher altitude 
regions are always filled with dense pair plasma. 

Screening of the accelerating electric field during each burst of 
pair formation proceeds slightly different than in cases described 
in previous sections. Particles are injected throughout the whole 
length of the gap and not at the gap's one end. Simulations shown in 
Fig. 25 have two "centers" within the gap where screening starts. In 
Fig. 26 we show in more detail than in Fig. 25 how the accelerating 



electric field is being screened by newly created pair plasma about 
the time t = 9.816. Fluctuations of charge and particle number den- 
sities are significantly less pronounced than those in Fig. 19 be- 
cause the accelerating field was due to mismatch of the charge den- 
sity Ijm/c — Tj\< \j GS \. The region of screened electric field spreads 
from two centers towards the edges of the gap. In Fig. 25 we show 
spreading of the second zone with the center around x = 0.2; the 
first "center" was at x — 0.1, the field screening started there earlier 
and electric field is already poisoned here. This spreading occurs 
again in the form of an electrostatic wave which phase velocity is 
larger that c and amplitude decreases with time. 

The position where the gap appears depends on the imposed 
current density and on the variation of the GJ charge density - the 
higher j m and the larger the variation of t] G! , the closer to the NS 
the gap starts developing. The accelerating electric field in this gap 
is smaller than that in the case with the constant GJ charge density 
from §6.2 because the number density of positrons is high and so 
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Figure 25. Screening of the electric field and final stages of pair formation in the flow with linearly varying GJ charge density and the imposed current density 
j m = 2 j Gs . The same quantities are plotted as in Fig. 11. Note that the time intervals between the first two snapshots are two times smaller that those between 
the rest of the plots. 



not the whole charge density j m /c is contributing to the electric 
field. If this electric field is not enough to produce pairs, then the 
gap grows up to high altitude until most positrons leave the domain; 
then a vacuum-like gap, with the charge density of the order of 
jm/c, starts developing at the outer end of the domain as in the case 
of constant GJ charge density from §6.2. 

The reason why variation of GJ charge density does not af- 
fect space charge limited flow with sub-GJ current density is the 
presence of the trapped particle "cloud" component which can ad- 
just to any charge density variation (assuming the current density 
remains sub-GJ), as was mentioned in §5. For flow with anti-GJ 
current density the accelerating electric field appears in the region 
depleted of electrons, which on average move toward the NS. In 
this case regions closer to the NS always have the larger number 
density of particles of the same sign as r/ GS and will become charge 
starved last, and so the gap develops at higher altitudes. For flows 



with the super-GJ current density when the absolute value of the 
GJ charge density decreases with altitude, a < 0, fewer positrons 
are necessary at lower altitudes to poison the electric field than at 
higher altitudes. Positrons are moving toward the NS and the low 
altitude region becomes charge starved after all the others and the 
gap develops at higher altitudes as in the case of the flow with a 
constant GJ charge density. 

Therefore, the regions with super-GJ flow can have two qual- 
itatively different behaviors depending on the value of the imposed 
current density and the variation of the GJ charge density. The cas- 
cade either starts close to the NS and high altitude regions are al- 
ways filled with dense pair plasma, or the cascade starts at high 
altitude and between the bursts of pair formation, the high altitude 
regions are charge starved with vacuum-like accelerating electric 
field which could give rise to an outer magnetosphere cascade zone. 
In the first case the period between successive bursts of pair forma- 
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Figure 26. Screening of the electric field and formation of superluminal electrostatic wave in the flow with linearly varying GJ charge density and the imposed 
current density j m = 2 j G] c. The same quantities are plotted as in Fig. 19. 

tion should be of the order the the gap's fly-by time and can be as 
short as ~ r pc /c ~ fractions of microseconds; in the second case the 
repetition rate will be larger than the fly-by time of the strong field 
zone, ~ R,/c ~ tens of microseconds. From our toy model we can 
not make quantitative prediction about the values of the imposed 
current density and magnitudes of the GJ charge density variations 
which separate those behaviors. 



6.4 On stationary cascades 

In all models previously considered in the literature, the size of the 
accelerating gap was limited by the opacity of the magnetosphere 
to the gamma-rays: pairs were injected when the optical depth to 
pair creation reached unity and those pairs screened the electric 
field. In our simulations with the constant GJ charge density, the 
electric field in the gap monotonically increases with distance, and 
the gap grows with time. That electric field could be screened only 
by injection of particles with the charge sign opposite to the charge 
sign of primary particles, and those particles, when injected, were 
accelerated in the direction opposite to the direction of motion of 



the primary particles by a strong electric field. They start emitting 
pair producing gamma rays after traveling a short distance, inject- 
ing pairs and so screening the electric field in their way. That led to 
the motion of the gap as a whole or to the gap closure, if the other 
end of the gap moved with subrelativistic velocity. 

In the case of a flow with varying GJ charge density, described 
in §6.3, the initial field was not a monotonic function of the dis- 
tance, but the gap was growing and, again, the only way to limit 
the gap and screen the electric field was by injection of particles 
of the opposite sign. Once injected those particles ultimately de- 
stroyed the gap, since their flux exceeded that of the primary beam 
- stationary equilibria with counter streaming particles of the op- 
posite sign exist only when the counter streaming flux is less than 
the primary flux, as in the Arons & Scharlemann (1979) model. 
Stationary particle acceleration and pair creation would be possible 
if pairs are injected (mostly) outside the gap, so that there will be 
not too many particles with charge sign opposite to that of primary 
particles within the gap. 

For example, the stationary flow model of Arons & Scharle- 
mann (1979) constructs a finite length acceleration zone bounded 
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Figure 27. Plasma flow with quasi-stationary cascade. The imposed current density is j m = 1.059/gj. it is "fine-tuned" in such way that the accelerating 
electric gap near the NS surface does not disappear. All other physical parameters are the same as in the case shown in Figs. 20,21. The same quantities are 
plotted as in Fig. 11. Note that the time intervals between snapshots are not equal. 



by the stellar surface below and the abrupt onset of pair creation 
above, which screens the electric field in a thin layer (the pair for- 
mation front, PFF). The space charge limited (electron) beam ex- 
tracted from the surface provides the primary particles that radi- 
ate pair creating gamma rays. Almost all of the pair creation oc- 
curs above the PFF. A small amount of trapping within the tran- 
sitional PFF layer sends charges of sign opposite to the primary 
beam (positrons) back down toward the surface with flux small 
compared to that of the primary beam, while extra electrons are 
added to the beam, restoring the charge density to equality with 
the Goldreich- Julian density. This restoration causes the screening 
("poisoning") of the accelerating field. In that model, E — > is 
imposed as a boundary condition, using the conclusion that dense 
pairs («+ » I'/cj/el) flash into existence at a rather well defined 
height - a physically correct conclusion when curvature radiation 
dominates the gamma ray emission and one photon absorption in 



the magnetic field is the major opacity source, since each primary 
beam particle emits many pair creating gamma rays and opacity is 
a rapidly varying exponential function of the photon energy. The 
application of that boundary condition, along with the space charge 
limited flow condition, has the consequence that the charge density 
of the beam has a unique value ^bcam- Since j = C7; bcara , this model 
becomes discrepant with the magnetospheric current density j m , in 
general. 

We tried to find an Arons & Scharlemann like solution with 
the acceleration gap limited by the pair formation front. For the lin- 
early varying GJ charge density given by eq. (16) and model param- 
eters from §6.3 we explored the parameter space by running differ- 
ent simulations for different values of the imposed current density 
starting from j m = 2 and decreasing it up to j m =s when pair 
production ceased. In each subsequent simulation with smaller j m 
the position of the first pair injection moved further from the NS 
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because the accelerating electric field was also smaller. In all cases 
when pair injection occurred inside the gap the gap boundary starts 
moving - we never saw a stationary PFF (nor did we see a station- 
ary PFF in any other simulations). So, the conclusion implied by 
the simulations is that the time asymptotic state is not a steady flow 
similar to Arons & Scharlemann's. Even though acceleration in the 
charge starved gaps is limited by pair creation with thin boundaries 
between pairs and quasi- vacuum, reminiscent of Arons & Scharle- 
mann's PFFs, the gaps move, either up into the magnetosphere or 
down toward the star - there is no truly stationary state, and fully 
developed limit cycles appear to be the answer, at least in ID. 

We now exhibit a novel steady flow model, which takes ad- 
vantage of the properties of non-neutral beam flow when rj G1 /B is 
nonuniform, that does not require poisoning by the pairs to produce 
E = at the gap's top. This model exploits the variation of the 
charge density with distance, in the case when the current density 
is larger than the local GJ current density only up to some height 
h„ then at distances larger than h s - where the non-neutral flow 
become sub-GJ - the electric field will be screened on the scale of 
the local Debye length, as discussed in §5. In such a situation the 
accelerating electric field exists only up to the distance h„ above 
which it is shorted out by a mixture of trapped electrons - those 
with two turning points in their orbits - and free pairs. 

We assume the GJ charge density varies linearly with the dis- 
tance, eq. (16), as in §6.3. By hypothesis, the flow is steady and the 
current density is equal to the imposed current density j m = 
£ > l.(seeeq.(12)).If£-l < a, at some point the non-neutral beam 
flow becomes sub-GJ. From the Poisson equation for the electric 
field (6), in the absence of pairs we have 



E s = Amf QS L 



(f-D : 



a I x 
2\ 



At the distance 



(19) 



(20) 



E s changes sign as the flow become sub-GJ. Identifying h s with the 
gap height yields the gap's potential drop to be 



(f-D 3 



(21) 



3 — a<- 

increasing in proportion to the imposed current density. If the re- 
sulting potential drop in the gap is large enough for particles to 
inject pairs within the gap with number large enough to form a 
positron back flux larger than the primary electron flux, the flow 
will be non-stationary, like that described in §6.3. But if the po- 
tential drop is such that copious pairs will be injected at distances 
larger than h s , positrons will be subjected only to a fluctuating, rel- 
atively small electric field, which sustain the cloud component at 
altitudes higher than h„ where the flow is sub-GJ. Only a fraction 
of the positrons will be advected into the gap, and, if this fraction 
will be much smaller than the GJ charge density the flow in the gap 
will be not disturbed enough to make it non-stationary. 

As an example of a flow with E(h s ) = in the non-neutral 
current flow, pair creation and non-neutral trapping, we show snap- 
shots in Fig. 27 of the flow with the same parameters as in §6.3, 
except with a carefully chosen current density which was set large 
enough to allow pair formation within the domain (at x < x B ), but 
small enough so that the particle back flux does not destroy the gap. 
The imposed current density in this model is j m = 1.059 j GS . Above 
the altitude x = h s « 0. 11 L the flow is sub-GJ and clearly has beam- 
cloud structure qualitatively similar to that in sub-GJ flows from §5 
(see phase portraits for electrons). The charged cloud of trapped 



electrons screens the electric field at x > h s . Primary electrons, 
when accelerated up to high energies, emit gamma-rays which de- 
cay into electron positron pairs above the gap. 

Both secondary electrons and positrons mix with the cloud 
component and some positrons are advected toward the gap. The 
gap shrinks a bit and the resulting potential becomes smaller than a 
critical value necessary to sustain pair production. When the num- 
ber density of positrons drops, the gap gets bigger and the pair for- 
mation starts again. The gap, however, never disappears and the 
cascade is close to a stationary configuration. This configuration 
is quite sensitive to the imposed current density; large gap fluc- 
tuations appear at imposed current densities only a few per cent 
larger than j m = 1.059/",. Our model, however, has an exagger- 
ated charge density gradient and a very high voltage and we expect 
that the mixing of positrons and their advection for real pulsar pa- 
rameters would be less efficient and, hence, there might be a larger 
parameter range (i.e. an interval of the imposed current densities) 
where this gap plasma flow can sustain quasi-stationary cascades. 

Such stationary configurations require the imposed current 
density being larger that the local GJ current density at the NS 
surface but then becoming less than the local GJ current density 
at some distance from the NS. Variations of the GJ charge density 
(after accounting for the decreasing magnetic field) would be not 
larger that =s 15% (Muslimov & Tsygan 1992; Hibschman & Arons 
2001) and, hence, the imposed current densities for stationary cas- 
cades are within the range \ j m - j GJ \ < 0.15, j m / j GS > 1. Taking into 
account that the resulting gap should not be very large to prevent 
pair injection within the gap, that range is even smaller. 

Fig. 1 makes clear that in the force-free theory, the current 
density mostly does not fall in this narrow range, lending support to 
the conclusion that limit cycle and trapping behavior are the generic 
physics for polar flow, with and without pairs. 



7 DISCUSSION 

These calculations, and those reported in Timokhin (2010), were 
designed to investigate the theoretical issue of how the pulsar mag- 
netosphere with current flow determined by the magnetosphere's 
global dynamics couples to the neutron star through the polar cap 
beam acceleration and pair creation zone. This problem was in- 
vestigated 30-40 years ago, primarily using order of magnitude 
estimates and analytical, steady flow models of charge separated 
beams, the accelerator was modeled as having voltage fixed by pair 
creation and geometry. That led to a determination of the charge 
density in the accelerator with the current then simply being charge 
density times the speed of light 10 . No effort was made to show that 
the current density estimated actually matched that required by the 
global structure - the models yielded currents with the correct or- 
der of magnitude, and in the absence of actual solutions for the 
global structure of the magnetosphere, that answer was regarded as 
good enough, although skepticism was expressed that operating the 
accelerator model as having fixed voltage actually could produce 
the correct answer (Mestel (1999) and references therein). The cal- 
culations reported here show that indeed, when the accelerator is 
operated with current fixed rather than voltage, the polar cap accel- 
erator does work, and works in a fully time-dependent manner, as 



10 In differing degrees, the charge density was made up of counter stream- 
ing relativistic beams, leading to current density lying between j] al c and 

2//GJC 
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first conjectured by Goldreich & Julian (1969) and Sturrock (1971), 
but with behavior different from previously published models 1 1 

Our results show that when the acceleration and cascades are 
one dimensional, acceleration and pair creation exhibits limit cy- 
cle behavior, with cycle time somewhat larger than the relativistic 
flyby time over the length of the accelerator. When the ID approx- 
imation is realistic (very young pulsars), the length of the acceler- 
ator is less than the polar cap diameter r pc ss 144P~ 1/2 meters, sug- 
gesting quasi-periodicity in the acceleration and pair creation on 
time scales > 10~ 6 />~ 1/2 seconds. As was suggested long ago, such 
variations, if reflected in the radio emission, might be the origin 
of the microstructure in the radio emission (Ruderman & Suther- 
land 1975). If limit cycle behavior is robust with respect to multi- 
dimensional considerations, this dynamical behavior might turn out 
to be a useful model for microstructure (in contrast to non-linear 
optical phenomena in the transfer of radio waves through the pair 
plasma, for example). 

When charges are fully bound to the surface (solid surface 
with high binding energy and no atmosphere, Medin & Lai 2010), 
current flow always adjusts to the magnetospheric load through pair 
creation discharges, as was shown in Timokhin (2010). The lowest 
energy particles in the pair discharge exhibit transient trapping in 
the fluctuating electric field. Those reversals of particle momenta 
allow both current and charge densities to adjust to the presumed 
force-free conditions, for any value of the magnetospherically im- 
posed current density. 

At long periods/weak B, pair creation becomes impossible, 
the magnetosphere is a vacuum with no conduction current flow 
and spin down occurs through generation of strong vacuum waves 
(Pacini 1967). Pulsars enter this regime by crossing the pair cre- 
ation "death valley", theoretically expected for magnetospheric 
voltages O m = V^V^ ~ 10 13 (/>/10- 15 ) 1/2 i>- 3/2 V < 10 12 = 
cD death V (a restatement of the death line based on curvature 
gamma ray emission and one photon magnetic conversion esti- 
mated by Sturrock (1971)). That this death line describes the disap- 
pearance boundary of radio pulsars in the P, P diagram reasonably 
well (see e.g. Fig. 15.1 in Arons (2009)) provides a strong hint that 
low altitude, polar cap pair creation with high energy beam acceler- 
ation and one-photon magnetic conversion of the curvature y-rays 
has something to do with pulsar activity. 

A warm, quasi-neutral atmosphere plasma atmosphere over- 
lying the magnetized ocean and solid crust, with no restrictions 
on charged particle motion along B, is the likely surface state of 
most pulsars, either because of the residual heat of the neutron star, 
or because of polar cap heating by precipitating particles from the 
magnetosphere (from local pair discharges or from the return cur- 
rent flow). The upper atmosphere then can freely supply charge, in 
a manner very similar to a space charge limited current flow from 
the cathode of a classical vacuum tube. Our simulation results on 
space charge limited flow with (and without) pair creation reveal a 
variety of noteworthy new features. 

We have made the first determination of the polar accelera- 
tor's behavior under conditions where the current density, not the 
voltage, is held fixed. We found that space charge limited flow is 
not universally high energy. Emission of pair creating gamma rays, 



11 Following (Shibata 1991), when the pulsar is not near "death valley" 
in P - P space, the load inserted into the magnetospheric circuit by the 
pair creation discharges creates a small perturbation of the magnetospheric 
current j m , thus justifying our treatment of j m as an external parameter in 
the description of the discharges. 



does not occur on all polar field lines, even when O m > Odcath- The 
force-free model of the magnetosphere suggests the polar current 
flow includes sub-GJ flow regions with < j/ j G1 < 1 , where the 
low energy current flow (with the Lorentz factor y beam < 3) with 
no (curvature) gamma ray emission and no pair creation (inverse 
Compton gamma ray emission and y-y conversion to pairs are both 
negligible), with adjustment of the current density and charge den- 
sity to the force free values through formation of a trapped particle, 
non-neutral hanging charge cloud. This kind of quasi-stationary 
flow occupies most of the polar flux tube for the aligned and al- 
most aligned rotator, but progressively disappears as the obliquity 
increases. The force free model also requires regions of distributed 
return current flow, jj j al < 0, and, as the obliquity increases, re- 
gions of super-GJ current flow, j/ j a! > 1 , occupying most of the 
polar flux tube as the obliquity approaches 90°. Both the return 
current and super-GJ regimes exhibit unsteady, high energy cur- 
rent flow (an unsteady beam), driving discharges copiously creating 
pairs whose character is similar to that encountered when the at- 
mosphere is absent and the surface is a strongly bound solid. These 
regions, as projected onto the polar cap, are shown in Fig. 1. In the 
special case of j m / j° s being slightly larger than unity, the plasma 
flow can sustain quasi-stationary cascades, and a new class of such 
models was described in §6.4. But generally speaking, such station- 
ary regimes represent a singular case, rarely if ever achievable by 
magnetospheres described by the force-free model. 

The only previous quantitative study of time-dependent cas- 
cades in space charge limited flow is that of Levinson et al. (2005). 
They used a two-fluid model and concluded that chaotic pair for- 
mation takes place throughout the whole zone where pair formation 
is possible. Our simulations do not support their conclusions. Parti- 
cle trapping plays a significant role in adjusting of the plasma flow 
to the required current density and, hence, the plasma can not be 
adequately represented as two fluids (electrons and positrons) each 
with its own unidirectional velocity. This two fluid representations 
introduced additional rigidity into the system and this, in our opin- 
ion, led to formation of strong chaotic electric field everywhere in 
the domain. 

The most recent analytical treatment of the problem was pre- 
sented in Beloborodov (2008). Our results support Beloborodov's 
(2008) general conjecture about the character of plasma flow in the 
force-free regime, that the sub-GJ flow is low energy and pair for- 
mation is possible only in super-GJ or anti-GJ current flow. The 
simulations differ from his qualitative picture of what would hap- 
pen in several respects: a) the low energy flow with trapped parti- 
cles does not retain the spatial oscillations of the cold flow, replac- 
ing those by the two-component beam/cloud structure - the warm 
beam has non-oscillatory velocity, while the trapped particles, not 
included in his picture, are those with two turning points in their 
orbits; b) bursts of pair cascades repeat after longer time than h/c, 
h being the gap height; and c) in most cases the gap is not destroyed 
but moves as a whole, usually relativistically. 

The discharges were given a novel treatment. After the semi- 
nal efforts of Ruderman & Sutherland (1975) but prior to our work, 
models universally incorporated their assumption that when pair 
creation is copious (many convertible gamma rays per primary high 
energy particle, as is the case in the curvature radiation cascades 
typical of young, high voltage pulsars), at and above the height 
where pairs start appearing, E\\ drops to zero and stays that way 
at all greater altitudes. Arons & Scharlemann (1979) formalized 
this into the dynamics of a pair formation front (PFF), showing 
that when pair creation is copious, the transition between a charged 
starved region where E\\ + and the pair dominated region where 
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£|l can be set equal to zero is thin - that structure necessarily in- 
volves one sided particle trapping, under pulsar conditions, thus 
causing the formation of counter streaming components in the cur- 
rent flow. Our solutions exhibit this transition between the dense 
plasma and the quasi-vacuum regions outside. But, since E\\ was 
obtained from a dynamical field equation, we nowhere assumed 
E\\ had to be zero everywhere outside (above or below) a pair 
discharge. Instead, the pair creation itself creates the polarizable 
plasma that dynamically resets Ey to zero inside the discharge, 
while outside the electric field was allowed to float to whatever 
is dynamically required. The resulting models thus exhibit macro- 
scopic intermittency - a pair discharge shorts out Ey, then as the 
cloud of pair plasma moves out into the magnetosphere or toward 
the surface, a gap reforms in which residual particles accelerate, 
radiate gamma rays and form a new discharge. The fact that pair 
creation opacity declines with distance from the neutron star was 
modeled by setting the magnetic field equal to zero in the upper part 
(typically the upper 30%) of the simulation domain. That Ansatz 
forces the one photon pair creation opacity to be zero in the upper 
part of the domain. 

It must be said, however, that the resulting model of the in- 
termittency is very simplified - probably over-simplified. The rep- 
resentation of the "magnetosphere" as a region of negligible pair 
creation optical depth to gamma rays emitted from just above the 
polar cap is the right general idea, but the ID character of our sim- 
ulations almost certainly leads to an overemphasis on the coher- 
ence of the intermittent discharges - they form coherent slabs of 
pair plasma, separated by coherent slabs of quasi-vacuum. In re- 
ality, the low opacity region where discharges must end appears 
at heights ~ R t » r pc = polar cap and polar flux tube transverse 
radius. Then the acceleration zone left behind a discharge created 
plasma cloud is long and skinny, with the narrow "wave-guide" 
geometry (and the physics of the polar flux tube's boundary) play- 
ing an essential role in the nature of the quasi-vacuum field. Under 
these circumstances, it is unclear and unknown whether the dis- 
charges would preserve anything like the coherence exhibited in 
the results reported here - formation of multiple "lightning bolts" 
is perhaps more likely, simultaneously coexisting within the polar 
flux tube, the electric fields of the discharges affecting the dynamics 
of their neighbors, causing the discharges to influence each other - 
a complex dynamical system quite likely to be chaotic . 

In addition, the extent of the gaps formed between discharges 
requires consideration of all the sources of pair conversion opacity, 
not only one photon conversion in the B field. Even when magnetic 
conversion drops to zero, y-y conversion with gamma rays collid- 
ing with soft photons from the atmosphere (either heated polar cap 
or overall warm surface, if the star is young enough) can provide 
discharge initiating pairs within quasi-vacuum gaps, limiting the 
extent of such gaps to less than what occurs when the opacity is 
set equal to zero. These issues require multi-dimensional model- 
ing of the accelerator and the photon transfer, as well extending the 
radiation transfer model, improvements required before the possi- 
ble direct consequences of low altitude discharges for observations 
can be addressed. Those consequences are the pair multiplicity of 
the outflow, the heating of the polar caps by the discharge compo- 
nents that move toward the neutron star, and possible direct col- 
lective emission of photons by the time dependent currents in the 
discharges. 

The multiplicity (k+ = n pa j r /n G j) within individual pair cre- 
ation bursts show the traditional results - when a discharge oc- 
curs, K ± S3 7beamra+C 2 /£ci,rvature ~ 10 3 (ficurvaturc = energy of pair 

producing curvature photons), somewhat enhanced by the conver- 



sion of the synchrotron photons in the cascade - thus, multiplicities 
within individual bursts up to ~ 10 4 are observed. Intermittency re- 
duces the overall multiplicity of the outflow - the magnitude of 
this reduction is hard to judge without a multi-dimensional model. 
The intermittency reduction is offset by continued pair creation as 
the plasma clouds and their high energy leading edges drift out to 
heights of order the stellar radius and more, emphasizing again the 
needs for a multi-D treatment. However, the multiplicity is unlikely 
to be as high as is inferred from pulsar wind nebulae (Bucciantini 
et al. 2011), so long as the simple, star centered dipole B field 
model is retained. More general B fields, with smaller radii of cur- 
vature of photon orbits with respect to B, can lead to enhanced pair 
creating opacity and pair yield - the offset dipole model (Arons 
1998; Harding & Muslimov 2011) is the most plausible magnetic 
anomaly model that can be consistent with the dipolar character 
of the polar flux tube revealed by radio observations (Rankin 1990; 
Kramer et al. 1998). Whether sufficiently large multiplicities can be 
obtained while respecting the geometric constraints obtained from 
the radio polarimetry is under investigation. 

Polar cap heating and consequent soft X-ray emission pro- 
vides another possible observable that can constrain the model. In 
common with all discharge models, roughly half the energy in each 
discharge is deposited in the crust below the atmosphere and ocean. 
If intermittency is neglected, existing observations provide strong 
challenges to all polar cap discharge models, including ours. In- 
termittency reduces the average energy flux. Whether that effect 
allows this polar discharge model to survive confrontation with 
observations (or, better, prove useful in modeling such observa- 
tions) remains to be studied. Simulating the particle back flux from 
a discharge is a particular challenge, requiring resolving the De- 
bye length in the pairs, which rapidly decreases as the pair density 
grows. 

Our results imply discharges occur on some of the polar field 
lines for all inclinations. These discharges incorporate time de- 
pendent, quasi-coherent currents on microsecond and shorter time 
scales. It has not escaped our attention that such fluctuations might 
be a direct source of radio emission from the low altitude polar flux 
tube, a region strongly suggested as the site of the radio emission 
by the radio astronomical phenomenology. Although the electric 
fields in our one-dimensional model are wholly electrostatic, there- 
fore cannot leave the plasma, in a more general multi-D setting 
which has substantial spatial inhomogeneity transverse to B, the 
field components parallel to B have accompanying components E ± 
perpendicular to B which make the fluctuations fully electromag- 
netic. Therefore, these electrostatic spectra may be representative 
of electromagnetic fluctuations which can leave the plasma, and 
the pulsar. From the point of view of simulations, a multidimen- 
sional, electromagnetic treatment is required in order to investigate 
this possibility 12 . 

In the frame of our model we can provide only estimates of the 
energy available for such directly excited waves. If the electrostatic 
oscillations could form an electromagnetic wave which leaves the 
plasma then the energy carried by the wave would be of the order 

W r ~ cnrl . (22) 



Fawley (1978) used a linear response analysis to study this possibility, 
in analytic model of discharges when strong surface binding suppresses fee 
emission. The work described in this paper is the first to study fully non- 
linear discharges in the free emission regime. Whether this idea for radio 
emission will lead to a useful model remains to be studied. 
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Figure 28. Estimated energy flux W in plasma waves for sub-GJ flows as a 
function off = j m / jci- W is normalized to Wm&BTiP 1 , 



Both kinds of flow, with or without pair formation, cause fluctuat- 
ing electric field. Fluctuating electric field in the sub-GJ flow turns 
out to be too low to provide the energy even for the radio emission. 
The electric field for the sub-GJ flow shown in Fig. 8 is normalized 
to Eq = |77gj|jT/Id,gj; in this normalization the estimates for W T is 
given by 



4.8- 10~ 9 < E 1 > W mi B\ l 2 P 2 



(23) 



where E is the normalized electric field and W m d is magnetodipolar 
energy losses 



W m 



BlR b t a. 4 

4c 3 



(24) 



In Fig. 28 we plot the estimated values of W T = W t /B~[lP 2 W m 4 as 
a function of £ = j m / j a . It is evident that even for the Crab pul- 
sar, known for its very low radio efficiency in terms of spin down 
energy losses ~ 10~ 8 , the energy in electrostatic oscillations of the 
cold flow can not power the radio emission. For plasma flows with 
pair formation, screening of the electric field proceeds similarly to 
cascades in the Ruderman-Sutherland model studied in Timokhin 
(2010). We observed formation of superluminal plasma waves dur- 
ing discharges in space charge limited flows, similar to what was 
found for strongly bound surfaces. The range of plasma oscilla- 
tion wavelengths was rather broad, similar to what was observed 
in discharges discussed in Timokhin (2010). The electric field for 
plasma flows with pair formation discussed in §6 is normalized to 
Eq = \Tia\nrp.; in this normalization the estimates for W t is given 
by 



1 



< W > W m 



(25) 



The amplitude of the (normalized) oscillating electric field during 
the screening phase of cascade development shown in Fig. 19, 26 
is ~ 0.1 - 0.01 and the energy in such oscillations is more than 
enough to power the radio emission. So, from the energetics point 
of view, it seems that the plasma flows with pair formation studied 
here are candidates for powering the pulsar radio emission through 
direct radiation by the discharge currents . 

Finally, our work may have implications for the formation of 
the particle accelerators in the outer magnetosphere required to ac- 
count for the gamma ray pulsars. Low voltage, sub-GJ current flow 
may have a Holloway (1973) "problem", in that the non-neutral 
cloud and low energy, non-neutral current-carrying beam cannot 



cross the null surface where CI • B = - in contrast, field lines pop- 
ulated with dense plasma from discharges have no difficulty with 
plasma flowing across the null surface, the pair plasma clouds are 
quasi-neutral and easily adjust to the locally required charge den- 
sity. Field lines passing from the low altitude trapped cloud to the 
outer magnetosphere that go through the null surface might form a 
physically self-consistent gap reminiscent of the earliest proposals 
for "outer gaps", with potential for outer magnetosphere discharges 
that could provide a model for the observed gamma ray emission. 
This issue also requires multidimensional modeling. The caustic 
formation exhibited by models for gamma ray pulses suggest that 
if such gaps can form, they are localized (by pair creation?) to re- 
gions close to the flux bundles where return currents flow. 

Another possibility for particle accelerators in the outer mag- 
netosphere exists in the regions carrying the return current. At some 
time between the bursts of pair formation in the polar cap the elec- 
tric field at large latitudes can not be screened by pairs created in 
the polar cap. This field can accelerate particles and give rise to 
pair creation via y-y process. As we pointed out at the end of §6.1, 
an outer magnetosphere cascade zone might exist along the mag- 
netic field lines carrying the return current. This might be a very 
intriguing possibility in view of the observational evidences that 
the spectrum of pulsar gamma-ray emission does not have a super- 
exponential cut-off, which should be present if there were absorp- 
tion of gamma-rays in the magnetic field. Moreover, modeling of 
pulsar gamma ray light curves suggest that the gamma-ray emission 
originates from the regions close to the boundary between the open 
and closed magnetic filed lines (e.g. Bai & Spitkovsky 2010). The 
return current regions for a broad range of pulsar inclination angles 
are close to that boundary or/and are within it, e.g. in Fig. 1 the 
return current is flowing in and around the current sheet for pulsar 
inclination angles a < 30° 



8 CONCLUSIONS 

Our principal conclusion is simple - pair creation can occur at 
pulsar polar caps, but (almost) always in the form of fully time 
dependent current flow (microsecond time scales for the variabil- 
ity). That time dependence with pair creation allows the current 
to adjust to any magnetospheric load, while simultaneously allow- 
ing the charge density to adjust to the requirements of the force 
free magnetosphere. We have also shown that a substantial fraction 
of the open field lines (fraction decreasing with increasing obliq- 
uity) solve the current flow problem with a low energy, non neutral 
beam carrying the current co-existing with a non-neutral, electri- 
cally trapped particle cloud. This is an essentially time independent 
local solution. Exploring the consequences of these new results for 
global theory and observations requires extending the calculations 
to multidimensional, electromagnetic models for the accelerating 
electric field, and perhaps to background magnetic field models 
more general that the star centered dipole geometry used here in 
the choice of the magnetic radius of curvature that enters into the 
pair conversion opacity. 
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APPENDIX A: ONE-DIMENSIONAL 
ELECTRODYNAMICS OF THE POLAR CAP 

Let us decompose the magnetic field in the acceleration zone into 
two terms B = B + SB, where B is the global time-averaged 
field (the "external" field for our local problem, = -B magnelosphere ) - 
and SB is the fluctuating magnetic field due to local currents in the 
system caused by charge motion in the acceleration zone. 

In our ID approximation the characteristic size of the accel- 
eration zone in longitudinal direction, along magnetic field lines, / 
is much smaller than its characteristic width w (the latter being of 
the order of r pc ). The characteristic timescale of electromagnetic 
field variations due to relativistic charge motion is t ~ l/c, as 
charges move along magnetic field lines. The characteristic scale 
of the global magnetic field variation in longitudinal direction Z mag 
is much larger than both / and w. The relation between these scales 

is / <<C W <K /mag. 

Electromagnetic field can be determined from Ampere's and 
Faraday's laws 



V x B + V x 6B 



An . IdE 
~ J+ c~dt 



Vx£ = 



(Al) 
(A2) 



1 dSB 

~~c~dF ' 

Combining Ampere's (Al) and Faraday's (A2) laws by eliminating 
the electric field we get an equation for the magnetic field SB: 



1 d 2 SB 



V 2 SB 



471 



-Wxj + W 2 B 



(A3) 



c 2 dt 2 

Only the perpendicular components of the magnetic field B ± ,SB ± 
affect the accelerating electric field E\\. Now using the perpendicu- 
lar component of eq. (A3) we estimate SB ± . 

The estimates of each of the terms in eq. (A3) are as follows. 
For terms with SB ± we have 



V Al p w 2 



l 2 



and 
1 d 2 SB 

3 dt 2 



SB ± SB ± 
(rr) 5 ~ 



(A4) 



(A5) 



The perpendicular component of the global magnetic field changes 
in the longitudinal direction on the scale / mag , in the lateral direction 
it changes on the scale w and, hence, 

(**o) ± ~ ^ + % ~ % • (A6) 

Current flows along magnetic field lines and its perpendicular com- 
ponent is of the order of (Bq i± /Bo) j ~ (l/p c ) j, and as the radius of 
curvature of magnetic field lines p c is much larger than any of the 
scales in our problem we have 

w I w p c I w 
Combining estimates (A4)-(A7) from eq. (A3) we have 



(Vxy\ 



(A7) 



SB ± 

w 



\wj \ c w I 



(A8) 
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The order-of-magnitude version of Ampere's law (Al) is 



Bo.x 8B L An g|| 

W W C ^ CT 



(A9) 



and from eq. (A8) it follows that the fluctuating magnetic field due 
to charge motion introduces only a second order term in l/w and it 
can be neglected. So, the Ampere law in ID has the form 



dt 



An 



(Vxfio 



ll) = -4*0'- 



7m). 



(A10) 



The same expression for the accelerating electric field through 
the current density as in eq. (A10) can be obtained from the Gauss 
law and charge conservation (see e.g. Timokhin (2010), Appendix 
A), in that case j m emerges as an integration constant correspond- 
ing to the time-average current density flowing trough the system. 
In ID the system is essentially electrostatic, charges create only 
electric field, and naturally both the Ampere's law and the Gauss' 
law reduce to the same equation. The magnetic field perpendicu- 
lar to the background stellar B field is negligible, both due to the 
background magnetospheric current and due to the rapidly vari- 
able currents within the acceleration zone considered in this study 
- in particular, to lowest order in (l/p c ) 2 , the particle orbits are de- 
termined by the stellar field B and consideration of the full elec- 
tromagnetism of the discharges is not required to characterize the 
discharge dynamics. 



APPENDIX B: STATIONARY ONE-DIMENSIONAL 
SPACE CHARGE LIMITED FLOW 

Bl General equations 

Let us consider a stationary one-dimensional beam of equal par- 
ticles with mass m and charge q having the same sign as the GJ 
charge density q GS flowing from the surface of the NS at x = into 
the magnetosphere at x > 0. The current density of the beam is a 
fraction £ of the GJ charge density, 



j = Z joj ■ 

For the energy of particles in the beam we have 
mc 2 y + q<p = mc 2 yo + q<po , 



(Bl) 



(B2) 



where y = (1 - v 2 /c 2 )~ 1/2 is particle's Lorentz factor, v is particle's 
velocity, if> is the electric potential; quantities with the subscript 
refer to their values at NS surface, at x = 0. The electric potential c/> 
is given by the Gauss law 



dx 2 



-An(q - 7] 0] ) . 



(B3) 



where q = j/v is the charge density of the beam. We consider the 
case when v, y, q, and <p are functions of distance x, while q GJ is 
constant. 

Expressing c/> through y from eq. (B2) and q through j = vq 
we get an equation for the beam's Lorentz factor 



d 2 y _ Anq Gl q 
dx 2 mc 2 



7 



1 



(B4) 



The plasma frequency for a mildly relativistic plasma with the GJ 
charge density consisting of particles with charge q and mass m is 



(Anq Gl q\ 



1/2 



Up.GJ - 



\ m 



1 



(B5) 



and the Debye length (also the skin depth, since the characteristic 
velocity is c) of such a plasma is 

^d.gj = — • (B6) 

Introducing normalized distance s = x/A D as eq. (B4) become 

d 2 y 



^v^r 1 - (B7) 

The first integral of this equation is easily obtained by multiplying 
both part by dy/ds 

</;.' ' 
ds 



= 2[f(V^-V^)-<r^ + (£) o - ^ 

Under the conventional assumptions the space charge limited flow 
starts at the surface with zero velocity and the electric field is com- 
pletely screened, so y = 1 and (dy/ds) = 0; with these assump- 
tions eq. (B8) become 

\2 

(B9) 



It is more convenient to analyze the properties of the flow in 
terms of the spatial component of the 4-velocity (the normalized 
momentum p = yv/c). In terms of p eq. (B9) is 

(i)'. a i±l (ltfr .^n). 

The right hand site of eq. (B10) is positive if either £ > 1 or < 

£ < 1 and p < /? raax , where 

2? 



Pmdx — 



W 2 



(Bll) 



So, for f > 1 flow will accelerate monotonically, while for < £ < 
1 the momentum will oscillate in the range [0, Pmail- 

According to eq. (B10) dp Ids does not depend explicitly on 
the distance s and so its absolute value is the same for the same 
value of p; therefore for < f < 1 the function p(s) is symmetric 
around its maxima and minima and is periodic. 

B2 Ultra-relativistic flow 

For p » 1 we can neglect terms of the order 0(1 lp) and higher. 
Then eq. (B10) takes the form 



2[l+p(f-l)] , 



(B12) 



the solution of eq. (B12) is (cf. eq. B3 in Fawley et al. 1977) 



V2i • 



(B13) 



If £ ^ 1 the flow is continuously accelerating and its momen- 
tum at large s will grow as 



(B14) 



For f < 1 the flow is oscillatory and its momentum period- 
ically reaches p max (see eq. (Bll)). If 1 - f <sc 1 and p max » 1 
the flow is relativistic almost everywhere. The wavelength of such 
spatial oscillations is twice the distance between the points where 
p = and p reaches its maximum value. From eq. (B13) the dis- 
tance where p reaches it maximum is V2(l -^T 1 and so the period 
of spatial oscillation is (Beloborodov 2008) 

.v = 2V2(l-^r'. (B15) 
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Figure CI. Numerical implementation of boundary conditions for space 
charge limited flow. Injection of one numerical particle is shown, particle's 
center is marked by a cross. See text for explanations. 



B3 Non-relativistic flow 

Stationary space charge limited flow is non-relativistic at the begin- 
ning, close to the NS surface; it can remain non-relativistic if f <c 1 
(seeeq. (Bll)). 

For p <sc 1 we can neglect terms of the order 0(p 3 ) and higher 
and eq. (BIO) becomes the Cycloid equation 

v 2 

v — n 

(B16) 



(dp) 
ids. 



2^-p 



The solution of this equation in parametric form is (Beloborodov 
2008) 



p = £(1 - cosw p , G ,0 

s = £(«p, G ,f- sinw p , G ,0, 



(B17) 
(B18) 



where the time t is measured in seconds. The period of spatial os- 
cillations is then 



(B19) 



For small if the flow is always non-relativistic and eqs. (B17,B18) 
describe it accurately everywhere. For large f eqs. (B17,B18) are 
good approximation near the starting point of the flow for all £ and 
for £ < 1 also near periodically repeating stagnation points, with 
small values of p. 

For small s, near the flow's starting point, the time t can be 
eliminated from (B17, B18) and we get (cf. eq. (13) inFawley et al. 
1977) 



p =s 



1/3 



f 1/3,2/3 



(B20) 



APPENDIX C: BOUNDARY CONDITIONS FOR 
MODELING OF SPACE CHARGE LIMITED FLOW 

For modeling of the space charge limited flow in the polar cap one 
needs to make an adequate numerical model for an infinitely large 
pool of particles available at the NS surface in order to correctly 
simulate the space charge limitation condition. This is not a trivial 
task. We found the following procedure works well. 

The calculation domain of the length L is divided in M x equal 



numerical cells (a typical value of M x in our simulations is ~ a few 
thousand). The electric field E and current j are set at cell bound- 
aries i = 0,..., M x . For calculation of the current density we use a 
ID version of the charge conservative algorithm proposed by Vil- 
lasenor & Buneman (1992), when charged particles are represented 
by uniformly charged sheets with the width equal to the cell size Ax 
and the position of a particle is the position of the sheet's center. 
The fraction of the sheet passed through the cell boundary i dur- 
ing a time step determines the contribution of the particle into the 
current j, at that point. 

At each end of the calculation domain we have one "ghost" 
cell. The outside boundaries of the ghost shells are "ghost" points 
with indexes i = -1 and i = M x + 1. The equation (12) for the 
electric field is solved for points . . . M x . Particles can move into 
the ghost cells but when their positions are outside of the domain 
[-Ax/2; L + Ax/2] and they are moving outwards they do not con- 
tribute into the current density in the domain anymore and such 
particles are deleted at the end of each time step. 

The electric field at the ghost points is set to zero. We solve 
a ID initial value problem - the electric field in the domain is cal- 
culated from the values of the electric field at the same points at 
the previous time step and is not coupled to the electric field at the 
"ghost" points. The electric field inside the ghost cells, at the par- 
ticles' position, are obtained by quadratic interpolation using val- 
ues E , Ei (and £m v -i> E Mx , E Mx+i ) and, therefore, setting E 
at the ghost point to zero {E-\ = E Mx +\ = 0) smoothly reduces 
the electric field inside "ghost" cells toward their outer ends. Set- 
ting the electric field at ghost points at each time step to the values 
obtained as extrapolation of electric field values near the domain 
boundaries (e.g. using quadratic extrapolation from points 0, 1,2 
and M x - 2, M x — l,M x ) or to some non-zero values resulted only 
in higher numerical noise and did not change the system behavior. 

At the beginning of each time step we inject certain amount of 
electrons A' in j and equal number of heavy positive charged particles 
("ions") in the first "ghost" cell of our computation domain at the 
position slightly outside the center of the fist cell x m j = -(Ax/2+5), 
see Fig. CI. The momentum of each injected particle is sampled 
from a uniform distribution in the interval [-pmj,Pmj]. In this way 
we can model finite temperature of the particles on one hand and 
populate the domain with particles more uniformly on the other 
hand, as each injected particle after the first time step will have 
a slightly different position. If during this time step the electric 
field inside the "ghost" cell cannot move particles into the inter- 
val [-Ax/2; L + Ax/2], they do not contribute to the current density 
and will be deleted at the end of the time step. Depending on the 
value of the electric field either positive of negative particles will 
be "sucked" into the domain. 

We experimented with different values for N m j and found 
that usually after N- m j exceeds the critical value = 
2(j m / Qc)(cAt / Ax) by ~ 20 - 30% computational results stop de- 
pending on A'inj ; further increase of N mj results in higher numerical 
noise. iVJjJ, is twice the number density of particles with relativistic 
velocities which must be injected at every time step in order to pro- 
vide the required current density j m ; Q is the charge of a numerical 
particle, c is the speed of light. The factor 2 accounts for injected 
particles having negative initial momentum - most of them do not 
reach computational domain and are deleted. Having some injected 
particles with negative momenta results in slightly lower numerical 
noise as well as more realistically represents the finite temperature 
of the NS atmosphere (the equivalent of the warm cathode in the 
analogous vacuum tube and high current beam technologies). The 
computational overhead caused by such particles is negligible, as 
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yVjnj is orders of magnitude less that the total number of particles. 
We also experimented with different values for the time step and 
found that values of At such that Ax/cAt ~ 5 results in relatively 
low level of numerical noise due to discrete events of particle injec- 
tion; smaller At lead to larger numerical overhead, as the stability 
of the leapfrog scheme requires only that At < 0.5Ax/c. 



© 2012 RAS, MNRAS 000, 1-35 



